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Abstract 



The transport of heat that results from turbulence is a major factor lim- 
iting the temperature gradient, and thus the performance, of fusion 
devices. We use nonlinear simulations to show that a toroidal equi- 
librium scale sheared flow can completely suppress the turbulence 
across a wide range of flow gradient and temperature gradient val- 
ues. We demonstrate the existence of a bifurcation across this range 
whereby the plasma may transition from a low flow gradient and tem- 
perature gradient state to a higher flow gradient and temperature gra- 
dient state. We show further that the maximum temperature gradient 
that can be reached by such a transition is limited by the existence, 
at high flow gradient, of subcritical turbulence driven by the parallel 
velocity gradient (PVG). We use linear simulations and analytic cal- 
culations to examine the properties of the transiently growing modes 
which give rise to this subcritical turbulence, and conclude that there 
may be a critical value of the ratio of the PVG to the suppressing per- 
pendicular gradient of the velocity (in a tokamak this ratio is equal to 
qje where q is the magnetic safety factor and e the inverse aspect ra- 
tio) below which the PVG is unable to drive subcritical turbulence. In 
light of this, we use nonlinear simulations to calculate, as a function of 
three parameters (the perpendicular flow shear, qje and the tempera- 
ture gradient), the surface within that parameter space which divides 
the regions where turbulence can and cannot be sustained: the zero- 
turbulence manifold. We are unable to conclude that there is in fact a 
critical value oiqje below which PVG-driven turbulence is eliminated. 
Nevertheless, we demonstrate that at low values of qje, the maximum 
critical temperature gradient that can be reached without generating 
turbulence (and thus, we infer, the maximum temperature gradient 
that could be reached in the transport bifurcation) is dramatically in- 
creased. Thus, we anticipate that a fusion device for which, across a 
significant portion of the minor radius, the magnetic shear is low, the 
ratio qje is low and the toroidal flow shear is strong, will achieve high 
levels of energy confinement and thus high performance. 
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Chapter 1 

Fusion, Turbulence and Sheared Flows 



1.1 Fusion Devices 

If one could persuade a nucleus of deuterium to collide with a nucleus of tri- 
tium at sufficient velocity, they would fuse together to produce helium, a neutron 
and 17.6 MeV of energy Of this energy 3.5 MeV would be given to the helium 
nucleus, and 14.1 MeV would be given to the neutron | ij]. 

If one could convince around 10^^ such pairs of nuclei to collide per second 
(that is, a total weight of just 10 milligrams of fuel per second), one could generate 
a gigawatt of power§ A mere 100 tonnes per year of fuel would power the United 
Kingdom,^ and 3500 tonnes per year the whole worldly 

The only product of this reaction is helium, a harmless and useful gas. The 
reaction produces neither carbon dioxide nor any of the radioactive byproducts 
of a nuclear fission power station. Although a small proportion of the vessel 
surrounding the reaction itself would become radioactive (at a very low level 
compared to the products of a fission reaction) over the life of the power station, a 
careful choice of design would ensure that the radioactivity of all but a negligible 
portion of that material would decline to safe levels within 50 years of the power 
station being closed |5|. Of the two fuels required, deuterium is abundant in 
seawater, and tritium can be bred from lithium (making use of the neutron) which 
is also abundant in seawater |^. 

This, then, is the case for fusion: a method of generating power which would 
require a mere 3.5 kilotonnes of fuel to power the entire world for a year, which 
produces no carbon dioxide and no radioactive waste. The fact that this source of 

^ Assuming a reactor efficiency of 30% 

^Based on gross UK energy consumption in 2010 of around 220 million tonnes of oil equiva- 
lent (Mtoe), with 1 Mtoe equivalent to approximately 42 GJ [3]. 
^ Based on world energy consumption of 8400 Mtoe in 2008 jj]- 
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energy still remains, after more than sixty years of research, beyond the reach of 
humankind, is owing to the problem of confinement: the problem of containing 
the fuel. 

In order to achieve the rate of 10^^ reactions per second needed to generate a 
gigawatt of power, 16 milligrams of fuel would need to be heated up (in a reaction 
volume of 10 cubic metres) to 7 million °C@ At such temperatures, no material 
known would be sufficient to confine the fuel directly. Therefore, another way of 
confining the fuel — which at such temperatures becomes a plasma — must be 
found. 

The most promising method to date is the use of a magnetic field. The fuel 
particles in a fusion plasma are fully ionised; as a result of their charge the electric 
and magnetic fields exert a force on them according to the Lorentz force law§ 

F = Ze(E+-vxB], (1.1) 



c 

where Ze is the charge of the particle, E is the electric field, c is the speed of light, 
V is the velocity of the particle and B is the magnetic field. A single particle in 
a strong magnetic field with no electric field is constrained by this law to move 
helically along a magnetic field line (see Fig. ll.lT a)). It cannot move across it. 
This means that if one could design a set of magnetic field lines that formed a 
closed surface, the field lines would act like the bars of the cage — the particle 
would never be able to cross that surface and escape. A well-known theorem ^ 
states that the only practicable closed surface on which a set of magnetic field 
lines can lie (without ever crossing the surface) is a torus. Thus, a magnetically 
confined fusion device principally comprises a strong equilibrium magnetic field 
whose field lines lie on a set of nested toroidal surfaces called flux surfaces, each 
of which forms a cage barring the way of a hot particle attempting to escape con- 
finement (Fig. I1.2[) . There are two principle types of such devices: tokamaks and 
stellarators. Tokamaks are devices which are roughly axisymmetric around the 
centre of the torus, in which the toroidal magnetic field (the component of the 
field around the major circle of the torus) is generated by external coils, and the 
poloidal field (the component of the field around the minor circle of the torus) is 
generated by a current which flows toroidally about the plasma |8, §1. Stellara- 
tors are non-axisymmetric devices whose equilibrium field is generated almost 



entirely by external coils |ld,lll[] . This work focuses exclusively upon tokamaks. 

Assuming a reaction cross-section oflO^ m^s" and a core density of lO'^hii^^ 
^ Except when stated otherwise, Gaussian units are used throughout this work. 
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Figure 1.1: (a) The helical path taken by a charged particle along a straight mag- 
netic field line, (b) The path of a particle moving across magnetic field lines via 
collisions with other particles at the marked points (an illustration only since 
head-on collisions as shown would be vanishtngly rare; see Footnote [S]). 




Figure 1.2: The toroidal cage of magnetic field lines which confines the plasma. 



The absence of a working fusion device based upon the principle of a mag- 
netic cage is clear evidence that this cage is far from perfect. Two effects which 
cause a lessening of the confinement provided by the magnetic field are the effect 
of collisions and the effect of drifts. In a typical fusion device there is not one par- 
ticle in the vessel, but of the order of 10^^. These particles collide with each other 
frequently: about 4 x 10^^ such collisions^ occur a second per cubic centimetre in 
the scenario described above ( see equation (|9.12[) ). A collision between two par- 
ticles can cause them to cross the magnetic field lines, and so in this way, through 
multiple collisions, a particle can eventually escape (see Fig. ll.lT b)). This mode 
of transport of particles (and with them heat and energy) across magnetic field 

*A collision here means a sequence of interactions which will change the momentum of the 
particle by order unity. 
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lines is known as classical transport. 

Classical transport itself is a minor problem. The effect of collisions is made 
far worse by the existence of drifts. Drifts, which can result from the electric field 
(the E X B drift), and the shape of the magnetic field (the gradient and curvature 
drifts), cause a particle to spiral slowly across magnetic field lines |12|1 . Drifts by 
themselves are not fatal. A fusion device can be designed so that on average a 
particle which drifts out somewhat will eventually drift back |13tl . It is the com- 
bination of drifts and collisions that is deleterious: a particle can drift some way 
out, collide with another particle and enter a new orbit which allows that parti- 
cle to drift further. Eventually the particle escapes. The transport which results 
from the combination of drifts and collisions is known as neoclassical transport. 
Neoclassical transport has been extensively studied |14l] . It is possible to directly 
calculate for many simpler fusion devices a good estimate of how large neoclas- 
sical transport is [15]; for the case of tokamaks, it is found that, although much 
larger than classical transport, neoclassical transport is not fatal: even including 
its effects, a fusion device only approximately 1 m in minor radius would be vi- 
able Q. 

The reason for the failure to produce a working fusion device lies in an alto- 
gether different phenomenon: turbulence. 



1.2 Turbulence 

Turbulence, in essence, is formed of fluctuations: fluctuations in the density 
and temperature; fluctuations in the electric and magnetic fields. Fluctuations 
in the magnetic field, by bending the bars of the magnetic cage, can cause very 
large degradations in the confinement of the plasma; in the core of many fusion 
devices, however, they are fortunately kept to a minimum]^ Fluctuations in elec- 
tric field can lead to eddies which lie across the magnetic field lines and which 
can convect heat, momentum and particles out of the plasma. The rate at which 
heat can escape via the turbulence is an order of magnitude greater than that 
due to neoclassical transport. It is the heat loss due to turbulence that has ne- 
cessitated the building of ever larger fusion devices, in order to have a greater 
distance from the core to the edge of the plasma (and thus a greater distance for 
the heat to travel to escape). The large sizes made necessary by turbulence create 
a raft of new challenges, particularly that of finding materials that can withstand 

^At least currently, see e.g. Ref. IitIi , for the simple reason that the particle pressure is in 
general low compared to the magnetic pressure: this may change as this ratio is increased ever 
higher in future devices. 
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the enormous heat loads, neutron fluxes and mechanical stresses that larger de- 
vices entail. It seems clear, therefore, that a vital task is to And ways of reducing 
or eliminating this turbulence. 

The turbulence in a fusion device is driven by the steep gradient in pressure 
between the core and the edge. Fluctuations in electric field are driven unstable 
by the pressure gradient and the resulting eddies grow to large enough ampli- 
tude to interact with each other, leading to turbulence. 

The pressure gradient can be split into the gradients of temperature and den- 
sity of the ions and electrons. The temperature gradients of both the ions and 
the electrons and the densit y g radient of the electrons can all drive turbulence 



19 



20l, I21L I22I, I23I, I24I, I25I, l26l, I271, 128(1 but the temperature gradient of the ions 



in particular is a strong driver of turbulence and eliminating the turbulent trans- 
port that results from the ion temperature gradient is the main focus of this work. 
The turbulent transport be eliminated both by preventing the growth of fluctu- 
ations, and by shearing apart the large eddies that result from the turbulence. 
Fortunately, there is a phenomenon which can achieve both: a sheared flow per- 
pendicular to the magnetic field lines. 



1.3 Sheared Flows 



A significant body of experimental work has shown that an equilibrium scale 
flow that runs perpendicular to the magnetic field lines which has a steep enough 
gradient across the magnetic flux surfaces can reduce, or even eliminate, turbu- 
lence 12911 . Such flows are often referred to as E x B flows as they can be thought 
of as a drift which results from a strong equilibrium radial electric field. They 



can be found both in the edge (where they are believed |30l. 



31 



32 



33(1 to add to 



the formation of the edge transport barriers responsible for the high confinement 
mode 341) ^rid in the core, where they can lead to a general reduction in turbu- 
lence (|35[l or to the formation of internal transport barriers, regions in the core 



m. In 



with strongly sheared flows and excellent confinement properties |36(, 3?!, 
this work we consider sheared flows within the core of a fusion device. 

Numerical work has confirmed the effectiveness of perpendicular sheared 
flows in reducing both the levels of turbulence and the associated heat loss 



41 



[t(42[,[43(,[4^[45(,[46(1. However, if there is a component of the flow paral- 
lel to the magnetic field, its gradient can drive another instability, in addition to 
the ion temperature gradient instability, known as the parallel velocity gradient 
(PVG) instability ([471]. As shown in Refs. (HHEl, the PVG instability can start 
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to drive turbulence at higher flow gradients. This matters because in general, 
the equilibrium flow within a fusion device is toroidal, and the angle between 
the flow and the magnetic field, and hence the ratio of the destabilising parallel 
flow shear and the stabilising perpendicular flow shear is fixed by the magnetic 
geometry, in fact by the ratio B^/Bg where B^ is the toroidal magnetic field and 
Bg is the poloidal magnetic field. This ratio is typically much greater than one. 
Nonetheless, Barnes et. al. |45ll showed that even at B^/Bg = 7.8, i.e., with the 
parallel velocity gradient nearly 8 times bigger than the perpendicular velocity 
gradient, it was still possible to quench turbulence completely, provided that the 
flow shear was large enough (and the ion temperature gradient was moderate). 



0.4 
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Figure 1.3: A cartoon of the effects of perpendicular flow shear, based on data 
from Ref . showing the growth rate 7 and the heat flux Q vs the perpendicular 
flow shear 7^;. The behaviour may be loosely divided into four regions (see text). 
(The normalising quantities are the ion thermal velocity Vthi, the tokamak major 
radius R and the gyro-Bohm estimate of the heat flux Qgs (see (|7.2)) ). 



The situation can be thought of as a three-way competition between the ion 
temperature gradient which drives turbulence, the perpendicular flow shear 
which suppresses turbulence, and the parallel flow shear which drives turbu- 
lence once again. Considering Fig. 2 of Ref. 11450, of which we provide a cartoon 



in Fig. 11.31 we may divide the interaction of the these effects, as a function of the 
gradient of the flow, into four regions. At zero flow gradient, an ion temperature 
gradient above critical drives strong turbulence, leading to the large heat fluxes 
at the left of the figure. In the first region, as the flow shear increases, the perpen- 
dicular flow shear reduces the levels of turbulence dramatically; in the second 
region the parallel velocity gradient causes a slight resurgence of the turbulence; 
in the third region the perpendicular flow gradient once again dominates, and 
can, at low temperature gradients, completely suppress the turbulence; at still 
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higher flow gradients, the fourth region, the PVG is dominant and the turbulence 
is rekindled. 

At lower flow shears, the strength of the turbulence is strongly affected by 
the growth rate of the underlying instability, which in turn is strongly affected 
by the flow shear (Fig. O after Fig. 1 of Ref. ||45D, see also Ref. |l48fl). At zero 
flow shear, the perturbations grow strongly for temperature gradients above the 
linear critical temperature gradient for instability. As the flow shear increases 
from zero, the growth rate drops precipitously owing to shear suppression in the 
first region, and then (at low temperature gradients) rises owing to the PVG in 
the second. In the third region the behaviour of the growth rate diverges from 
that of the turbulence; we find that for all temperature gradients the growth rate 
drops completely to zero. It remains zero even in the fourth region, where the 
turbulence starts growing again. 

There are two points of great interest in the behaviour of the linear growth 
rate. The first is that it drops to zero at all temperature gradients in the third 
region. This behaviour is explained by Newton et. al. |49(1 as follows. At low flow 
shears, where there is a finite magnetic shear or twist in the magnetic field lines, 
it is possible for the magnetic shear to exactly cancel out the flow shear, provided 
the perturbations move along the magnetic field lines at a speed proportional to 
the flow shear divided by the magnetic shear. As the flow shear increases, this 
speed becomes higher than the sound speed in the plasma and it is no longer 
possible for the perturbations to move quickly enough; they are sheared apart by 
the flow and die away. 

The second point of interest is that it is possible to have strong turbulence 
even when the growth rate of the instability is zero. This is possible because 
even though the perturbations eventually die away, they can grow transiently 
for long enough for their energy to be scattered by nonlinear action into other 
perturbations which are still growing (as illustrated in a simple matrix model by 
Baggett 1 5^). Turbulence which exists where there are no linear instabilities is 
known as subcritical turbulence and will be discussed extensively in this work. 



1.4 The Zero Magnetic Shear Limit 

In this work we consider the effects of flow shear when the magnetic shear 
is zero. When the magnetic shear is zero the speed with which a perturbation 
would need to move in order to grow without limit is infinite for all non-zero 
values of the flow shear. Thus, all perturbations grow only transiently before 
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decaying. We find that the first and second regions as described above disap- 
pear and that all turbulence is subcritical. We find also that in the absence of 
growing modes the turbulent transport levels for any given flow or temperature 
gradient are always lower than at the finite value of magnetic shear (s = 0.8) 
considered in Ref. [45|] . We find, in addition, that at zero magnetic shear condi- 
tions are favourable for a bifurcation from low to high flow gradients, and from 
low to high temperature gradients, of a sort that may lead to the formation of 
a transport barrier. This agrees with the experimental observation that internal 
transport barriers tend to form in regions of low magnetic shear |51.l52]. 

The rest of this work is organised as follows. In Part [III we discuss the way 
that turbulence and the effects of flow shear are modelled. In Part|lll]we calculate 
the effect of the flow shear on the turbulent transport using nonlinear gyrokinetic 
simulations. We demonstrate the existence firstly of a large range of values of the 
flow shear and the ion temperature gradient where the turbulence is completely 
suppressed, and secondly of a transport bifurcation from low to high flow gradi- 
ents and from low to high temperature gradients. We present a parametric model 
of the turbulent transport and use it to investigate further properties of the high- 
gradient, reduced-transport state reached by the bifurcation, as well as the values 
of the heat and momentum fluxes (which, assuming steady state, are the control 
parameters of a real fusion device) for which we expect such bifurcations to occur. 

In Part |IV] we investigate the phenomenon of subcritical turbulence in more 
detail. We start by using a linear model in simplified geometry to calculate the 
way the strength of the transiently growing modes is affected by the parameters 
of the system. We identify three key parameters, the ion temperature gradient, 
the flow shear and the ratio q/e, which strongly affect the strength of these tran- 
sient modes. We then use nonlinear simulations to map out, as a function of these 
three parameters, the zero-turbulence manifold, the dividing surface between the 
area of that parameter space where the turbulence can exist and the area where it 
cannot. If we cast the problem of improving confinement as that of increasing the 
ion temperature gradient (and hence the core fuel temperature) as much as pos- 
sible without generating turbulence, we find that the optimum strategy is firstly 
to minimise the ratio g/e (as might be expected), and then to increase the toroidal 
flow shear to an optimum value which we calculate. 

In the final section of this work. Part |Vl we present our conclusions and our 
suggestions for further work. 
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Part II 

Modelling the Turbulence: A 
Concise Description of GS2 
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Chapter 2 
Modelling the Turbulence 



2.1 Introduction 

The aim of modelling plasma turbulence in a fusion device is to determine, 
for a given set of plasma properties, for example a given temperature gradient, 
how much heat, momentum or how many particles escape through the plasma as 
a result of that turbulence. The task of modelling the turbulence can be divided 
into two parts: deriving an equation which describes the way the turbulence 
evolves in time, and solving that equation. In general, the equations that govern 
plasma turbulence cannot be solved analytically and must be solved numerically. 

In this work we use the gyrokinetic equation to describe the evolution of the 
turbulence. The gyrokinetic equation is derived from the Vlasov-Landau equa- 
tion, which describes the evolution of a general plasma. The gyrokinetic equation 
is specifically designed for modelling the particular phenomena of microturbu- 
lence, for example the turbulence that is generated by the ion temperature gra- 
dient in a tokamak. It uses various properties of the microturbulence to simplify 
the Vlasov-Landau equation and make it easier to solve, in particular reducing 
the number of dimensions in the problem from 6 to 5. These properties are, in 
summary, that it evolves on timescales that are much faster than the evolution of 
the plasma equilibrium, but much slower than the plasma frequency or the gy- 
ration of the particles about the magnetic field lines; that the turbulent structures 
are much smaller than the scale on which the equilibrium changes, and that the 
turbulent structures are strongly elongated along the magnetic field lines. The 
derivation of the gyrokinetic equation based on these and other assumptions has 
been extensively covered elsewhere [ jssll . and is only covered in summary here in 
Section IZ2l 

The tool which will be used to solve the gyrokinetic equation and calculate the 
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behaviour of the turbulence is the code GS2. We will use GS2 to calculate the be- 
haviour of the turbulence within a small region of the fusion device called a flux 
tube. To do this we make a set of further simplifications known as the local ap- 
proximation. In particular, we assume that the equilibrium density, temperature 
and flow gradients are constant across the flux tube. 



2.2 The Gyrokinetic Equation 

In order to model the turbulence, it is necessary to have an equation which 
describes how the turbulence evolves in time. In this work we will use the gy- 
rokinetic equation |5^ in the high-flow low-Mach limit. The derivation of this 
equation, including all the assumptions and approximations that go into it has 
been extensively covered, in particular by Abel et. al. |53l] and Sugama and Hor- 



ton |55|l . Here, we only provide a brief summary of the derivation, and refer our 



readers to Ref. 153(1 ^rid references therein for a comprehensive description. 

We consider a toroidal plasma with an axisymmetric equilibrium magnetic 
field B. The field lines lie on closed flux surfaces, which can be labelled by the 
poloidal magnetic flux tp contained within each surface (see Section |3] for a much 
more detailed discussion of the equilibrium magnetic field). The magnetic field 
can be written as 

B = V(j)xV^ + (2.1) 

where (p is the toroidal angle, R is the radial coordinate measured from the cen- 
tral axis of the torus (the major radius), and I/R = is the toroidal magnetic 
field. Owing to the fast motion of the particles along the magnetic field lines, 
equilibrium quantities such as the density and temperature are constant on each 
flux surface II 

We allow an equilibrium plasma flow u, of the same order as the ion thermal 
velocity 

vthi = \ —, (2.2) 

V m 

where Tj is the temperature of the ions and nii is their massj^ It can be shown 
that such a flow is toroidal, as any poloidal component will be quickly damped 



^ For the density, this is only true if the toroidal flow is subsonic, which is what we assume 
below, see Eq. 12.91 

^ The definition of the thermal velocity given in l|2.2b is chosen in accordance with the analyt- 
ical works [53, 15] cited in this thesis; maintaining correspondence with [45, 46], where the ion 
thermal velocity was defined as \JT\fmi, results in various factors of ^2 in the normalisations of 
the output from simulations. 
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|56l, 115(1 . The flow can then be expressed as u = uR^Wcj), where u is the angular 
velocity of the flow (which must be constant on a flux surface). 

The state of each species s of charged particles can be described by its distri- 
bution function, fs, which is, very loosely, the probability that there is a particle 
at a given location x and travelling at a given speed v. Thus, fs is a function with 
seven parameters; three spatial coordinates, three velocity coordinates, and time. 
The evolution of fs is given by the Vlasov-Landau equation: 

dfs , dv dfs . , 

where C, the collision operator, comprises the effects of collisions between parti- 
cles. The Vlasov-Landau equation could, in theory, be solved directly. However, 
since it is a six dimensional equation, involving scales which range from the tiny 
orbit size of particles gyrating around the magnetic field lines to the size of the 
device, its direct solution for a fusion device is impossible. In order to proceed, 
we make certain assumptions which lead to an equation which is in fact soluble: 
the gyrokinetic equation. 

Firstly, the Vlasov-Landau equation can be expanded by splitting fs into an 
equilibrium part Fs and a perturbation 5fs. 

fs = Fs + 6fs. (2.4) 

The latter is further split into averaged and fluctuating parts: 6fs = {6fs) + Sfs, 
so that 

fs = Fs + {6fs) + 6fs (2.5) 
where the angle brackets denote a spatial and temporal average over fluctuations 



|53|1 . Secondly, we assume that: 



the perturbations to the distribution function (and the background elec- 
tromagnetic fields) are small compared to the equilibrium values of those 
quantities; 

owing to the fact that particles can travel quickly along field lines but only 
slowly across them via collisions or drifts, the perturbations vary slowly 
along the field lines but quickly across them; 

the perturbations vary very quickly compared to the rate at which the equi- 
librium changes but slowly compared to the speed with which the particles 
gyrate around the magnetic field lines; 
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the spatial scale of the perturbations is small compared to the scale of vari- 
ation of equilibrium quantities. 



These assumptions give rise to the gyrokinetic ordering |54j, |53|1 which may be 
written: 



where k\\ and fc^ are the typical parallel and perpendicular wavenumbers of the 
fluctuations, 7 is the growth rate of the fluctuations, Vli and pi are the gyrofre- 
quency and Larmor radius of the ions, respectively, L is scale length of the vari- 
ation of Fs, and egk, the gyrokinetic expansion parameter, is small. Using the 
gyrokinetic ordering, it is then possible to show]^ that to lowest order, F, is a local 
Maxwellian in the frame of the equilibrium flow u: 

where Ss is the energy of the particle, nxg is the mass of species s, and Ug and 
are the equilibrium density and temperature of species s, respectively. It may be 
further shown that to first order in egk, 

6fs = — + hs {t, Rs, fis, Ss) , (2.8) 

^ s 

where ZgC is the charge of species s and (p is the perturbed electrostatic potential. 
The non-Boltzmann part of the distribution function, hg, the gyrocentre distri- 
bution, is independent of the gyrophase, that is, the exact position of a particle 
within its orbit around a magnetic field line. Thus hs is only a function of the cen- 
tre of that orbit Rg, the perpendicular velocity v± of the particle around that orbit 
(for which we use the coordinate fig = msv\/2B), the energy Ss of the particle. 



and the time, t 153(1 . The gyrocentre position Rg = r — b x v/Vtg, where r is the 
position coordinate, Vlg is the gyrofrequency of species s, the unit vector h = B/B 
is the direction of the equilibrium magnetic field, and B is its magnitude. Thus, 
the problem has been reduced from a six-dimensional to a five-dimensional one 
by only considering timescales longer than the fast orbit of the particles around 
the magnetic field. 

The turbulence can be affected both by the gradient of the flow and by its 
magnitude, the latter through the Coriolis and centrifugal drifts. Here we neglect 



^ss < 7 



With the additional assumption that the like-particle collision frequency Vgs satisfies e^^'-^ <C 
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the effects of the magnitude of the flow (these effects are described in detail in 
Refs. llsz, 58]) but allow the gradient of the flow to be of order the growth rate of 
the fluctuations, Rdco / dr ~ 7, as is necessary for the flow gradient to affect the 
turbulence non-negligibly |39|1 . This can be systematised by assuming that the 
magnitude of the flow is small but the gradient of the flow is large in a Mach- 
number ordering subsidiary to (|2.6|) , namely 

^ ^ . Rdu 1 ^ „, 

M<1, -— ~ — >1 (2.9) 

Vthi ^ dr M 

where egk < Af < 1. 

Under these combined assumptions, it can be shown that hg evolves according 
to the following gyrokinetic equation: 



— + u-V]hs+ [v\\b+VDs + ^bxV{ip)ji^ 



Ts 



c ? 



dFs msIv\\Fsduj 



bxV{ip)J,)■V^P, (2.10) 



where ()^^ is a gyroaverage (an average over a particle orbit) at constant Rs, 



V, 



Ds 



vffi-Vb+ -V\nB 



(2.11) 



is the magnetic drift velocity, v^^ = crA/(£s — fJ-sB) /rrig is the parallel velocity, a 
is the sign of f | , and c is the speed of light. We further simplify the problem by 
assuming purely electrostatic perturbations {5B = 0), so the system is closed by 
the quasineutrality condition: 



(2.12) 



where 0^, is a gyroaverage at constant r. 

In general, the gyrokinetic equation must be solved separately for each species 
s, and the resulting perturbed distribution function for each species hg entered 
into the quasineutrality condition to compute how each species affects the elec- 
tric field. However, in this work, we do not solve the gyrokinetic equation for 
electrons, but instead treat them using a modified Boltzmann response 1 55 

e{ip-lp) 



dMhe), 



(2.13) 



where the overbar denotes averaging over the flux surface. Thus, Eq. (|2.10|) is 
only solved for ions, s = i, and together with (|2.12[) comprises a closed system of 
equations for the evolution of h^. 
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Knowing hi, we can calculate the turbulent heat and angular momentum 
fluxes associated with a given minor radius and given values of the gradients. 



Qt = l^j d^v^rmv^^iki}^ (b X V<^) • Vr^ (2.14) 

= Q d^vm,{h,v)^ ■ V0i?2-| (h X V<^) ■ Vr^ (2.15) 
All that remains is to solve the equations (|2.10|) and (|2.12[) and find h^. 



2.3 An Overview of GS2 

In this work, the gyrokinetic equation derived above is solved numerically 
using the freely available simulation code GS2 |6^. GS2 is an initial value code 
which solves the nonlinear 5f continuum gyrokinetic equation in a small region 
of the plasma known as a flux tube, in which it is assumed that the gradients of 
equilibrium quantities are constant. It is one of the earliest in a class of codes 
which solve the 5f gyrokinetic equation in similar ways; other notable members 
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of this class include the codes GENE |27 

There is no one existing source which fully describes the exact equations and 
algorithms of GS2; as the majority of of the results contained within this thesis 
are generated using it. Chapters |3ll6] are devoted to describing how GS2 works J 
Here we give a brief overview of GS2 and the extensive, if fragmented, existing 
documentation. The main algorithms of GS2 may be roughly summarised as: 



the irnplicit algorithm which solves the linearised set of equations (|2.10[) and 
(I2J2H 



• the addition of the nonlinear term (the term proportional to 6 x V (v^)^^ ■ V/ij 
in (|2.10|) ) as a source in the linear equation, 

• the addition of collisions. 

The original exposition of the implicit algorithm that solves the linear equation 
in the absence of collisions is presented by Kotschenreuther et. al. [67], but, as 
is documented by Belli [6^, the actual implementation of the algorithm is dif- 
ferent. The clearest and most compact description of the linear algorithm as it 

• In the absence of electromagnetic perturbations. 

^ It should be noted that a simulation code like GS2 is continuously evolving. The description 
given in this work corresponds to revision 1555 of the trunk in the Subversion repository on the 
Sourceforge Gyrokinetics project page fsoh . 

^With the addition of Ampere's Law if magnetic perturbations are included. 
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Table 2.1: A summary of which parts of GS2 are documented in each source. The 
letters are defined as F — full: covered in full detail; C — concise: the whole func- 
tionality of the code is covered, but briefly; P — partial: only a (possibly large) 
part of the functionality of the code is documented; — old: the description does 
not correspond to the code in its current form. 



Source 


Linear Algorithm 


Nonlinear Terms 


Collisions 


Electromagnetic Terms 


Coordinates 


Geometry 


Boundary Conditions 


Velocity Space 


Flow Shear 


Normalisations 


Here 


F 


F 


C 




F 


P 


F 


C 


F 


P 


Beer et. al. [66] 










F 


P 


P 








Kotschenreuther et. al. [671 














P 


P 









Dorland & Kotschenreuther [68] 










P 


F 


P 






P 


Belli [6^ 


F 






F 














Barnes & Abel [70, 71J 






F 














P 


Numata et. al. [72] 


F 


C 


C 


F 








c 




P 


Barnes [73] 






F 






F 




F 




P 



currently stands is given in the paper documenting the AstroGK code |72[] . but 
as AstroGK is a straight-field-line simplification of GS2, it does not cover some 
of the complications that follow from the geometry, such as the parallel bound- 
ary conditions, which affect GS2. The inclusion of the nonlinear terms is noted 



in Ref. |27|] and described more fully in Ref . |72l] . while the collision operator is 



derived by Abel et. al. ||zQ] and implemented according to Barnes et. al. [71]. 
In addition to the main algorithms, there are additions and subtleties such as: 

• the definition of the coordinates, 

• the treatment of the geometry of the magnetic field, 

• the treatment of boundary conditions, 

• the layout of the velocity space grids, 

• the implementation of perpendicular flow shear, 

• the normalisation of input parameters and calculated quantities. 
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Again, the documentation of these is fragmented; the treatment of boundary con- 
ditions is discussed in Ref . |67l] , following the work of Beer et. al. [66]; the layout 
of the velocity space grids and the calculation of velocity integrals is described by 
Barnes IZst] and the implementation of flow shear is briefly described in Ref. [74]. 
Descriptions of coordinates, geometry and normalisations are giv en in an unpub- 
lished paper |6^ which was included as an appendix to Ref. |Z3]. A summary of 
the current documentation of GS2, as described below, is given in Table IZTl 

The rest of Part|II]is organised as follows. In Chapter |3j we describe how the 
geometry of the equilibrium magnetic field is dealt with, the choice of coordinates 
and the variables used within GS2 and the way they are normalised. In Chapter 
m, we describe the algorithms that are used to solve the gyrokinetic equation in 
the absence of collisions between particles. In Chapter|5j we describe the way that 
the velocity integrals of the distribution function (used to calculate, for example, 
the perturbed field generated by the particles) are treated, and how the effect of 
collisions between particles is included. In Chapter [6] we describe how the effect 
of a sheared equilibrium flow is incorporated within GS2. 
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Chapter 3 



Geometry 



3.1 Introduction 

Many of the terms in the gyrokinetic equation (|2.10|) depend on the shape or 



strength of the magnetic field; they are given in Table I3.1[ In a fusion device, 
the magnetic field can have a very complicated shape indeed, and the task of 
calculating terms that depend on it becomes equally complicated. Fortunately, in 
a well-behaved equilibrium in a fusion device, that is, in the absence of magnetic 
islands and stochasticity, the magnetic field lines have properties, namely that 
they are continuous and do not cross, that make them an excellent basis for an 
curvilinear coordinate system, with one unit vector along the magnetic field line, 
and two perpendicular to it. In any such coordinate system, many operators 
such as b ■ V can be expressed simply. Such a coordinate system also allows a 
simulation code to take advantage of the fact that structures in the turbulence 
tend to be elongated along a background magnetic field: this allows grid points 
along the magnetic field to be much more widely spaced. This saving in grid 
size is indispensable, and it is principally for this reason that field line following 
coordinate systems are used. 

The two tasks that must be faced, therefore, when dealing with the geometry 
are the definition of a set of field line following coordinates, and the calculation of 
those terms in the gyrokinetic equation which depend on the shape or strength 
of the magnetic field and which are not rendered trivial to calculate by such a 
coordinate system. The first task can be easily achieved (at least, once it has 



been done once: Ref. 11660); the second is Herculean. The reason for this is that 
while the choice of a field line following coordinate system renders most of the 
gyrokinetic equation easier to solve, the calculation of the remaining terms such 
as the magnetic drifts requires one to keep track, firstly, of which points in real 



3.2. Defining a Coordinate System 



Operator 


Description 


u ■ V 




Equilibrium Flow 


b- V 




Parallel Streaming 


6 X (6 ■ V6) 




Curvature Drift 


b X (Vln5) 




Gradient Drift 


(c/B) {b X V((^)^J 




Nonlinear Term 


(c/5) (6 X V(y.)^J 




Perturbed Flow (radial component) 






Background Gradients 






Perpendicular Wavenumber 


B 




The Magnetic Field Strength 



Table 3.1: The operators in the gyrokinetic equation which depend on the size or 
shape of the magnetic field, along with their physical significance or the physical 
significance of the terms in which they appear. 



space the grid points in this new coordinate system correspond to, and secondly, 
what the shape and strength of the magnetic field actually is at those points. 

3.2 Defining a Coordinate System 

When defining a coordinate system, we assume that the magnetic field vector, 
B, is known at every point. All coordinates will be defined with respect to this 
vector. We will now define our three basis vectors and our three coordinates. 
One of our basis vectors is, naturally, the magnetic field direction: b = B/B (we 
choose the corresponding coordinate below). The second follows from the fact 
that in any axisymmetric toroidal system, the magnetic field lines will all lie on 
nested surfaces, and each of these nested surfaces encloses a well-defined poloidal 
flux, that is, a magnetic flux going around the magnetic axis, which is defined as 
the integral of the poloidal magnetic field Bg = B- V6 across a surface of constant 
poloidal angle 6 (see Fig. 13.1)1 which extends from the magnetic axis to the flux 
surface itself: l] 

^ (r) = / dr' r'R{r)B{r) ■ Vd. (3.2) 



^ In fact, it is perfectly possible to have a non-axisymmetric system where the magnetic field 
lines still lie on well-defined closed flux surfaces (for example, a stellarator). In this case, we can 
use a more general definition of tp: 



2 \ 2 rV 



^ = ( — ) / dVB-Ve (3.1) 







where V is the volume enclosed by the flux surface. This formula reduces to (I3.2t for the case of 
an axisymmetric system. 
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The magnetic axis is at the centre of the nested flux surfaces, the point where 
Bg = and the field is purely toroidal. The minor radius r is the distance from 
the magnetic axis, and R is the major radius at r (see Fig.l3.1|). Since the quantity tjj 
is constant on a flux surface, and is therefore constant on any magnetic field line, 
it will serve as our second coordinate, and its gradient will serve as our second 
basis vector: i/j = VV'/ IV-?/']. 

For our third basis vector a, we require b ■ a. = (so that coordinate system 
follows the field lines) and b- {ax ijy) ^ {to prevent the coordinate system from 
diverging) at every point except at the magnetic axis. A coordinate a (with basis 
vector a. = Va/ |Va|) defined such that B = Va x VV^ and Va ■ Vip ^ B will 



satisfy these conditionso Kruskal and Kulsrud |76|1 showed that for a magnetic 



equilibrium composed of closed flux surfaces, we may write: 

a = (l) + q{^l))d + y{<P,d,ij), (3.3) 

where is the toroidal angle (see Fig. |3.1|l , is the poloidal angle, z/ is a func- 
tion which depends on the geometry and which is periodic in and 6, and the 
magnetic safety factor 



where il)t is the toroidal magnetic flux 



i 



i^t= £ dVB-V(j) (3.5) 

Lastly, for the coordinate along the magnetic field, GS2 follows Ref. [66 ] and 
many others in choosing the poloidal angle 9. This may seem counterintuitive: 
why not just choose a coordinate such as "the distance along the field line"? In 
fact, we note firstly that at fixed ip and a, any coordinate which is not fixed will 
serve as the parallel coordinate, and secondly that, at fixed many equilibrium 
quantities (such as B) can be written (in some limits) as quite simple functions of 
9. 



3.2.1 The Local Approximation 

In order to take full advantage of the separation between the scale of the tur- 
bulence and the scale of the equilibrium GS2 solves the gyrokinetic equation in 
a small region of the plasma known as a flux tube. A flux tube is designed to 

^ In this case, a is a surface potential (t^ . 

^ No equivalent of the Il3.2l l exists for the toroidal flux as in general the flux surfaces in a 
tokamak are not symmetric in 9. 
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Magnetic Axis 



R 



r 




a 



6 



Figure 3.1: Toroidal coordinates. Note that the signs of 9 and are not fixed (see 
Section |3.3[) , though the origin of 6 is fixed at the height of the magnetic axis on 
the outboard side. Note also that a^, the half radius of the flux surface at the 
height of the magnetic axis, is not in general equal to r, the minor radius (except 
in the very special case of zero Shafranov shift and circular flux surfaces). 

be large enough to include several turbulence decorrelation lengths in both the 
parallel and perpendicular directions (i.e. large enough that the turbulence in the 
centre cannot "see" the edge of the box) but small enough that high spatial reso- 
lution can still be achieved. Since the turbulent eddies themselves are elongated 
along the magnetic field line, a flux tube likewise follows the magnetic field lines 
and is elongated along them (see Fig. 1 of Ref. |66l]). Within the flux tube we 
define two local perpendicular coordinates, x and y, which measure the distance 
from the magnetic field line at the centre of that flux tube, located at (V^o, «o): 



X 



1 go ^ 
aBaPo 
1 dip 



{ip - ipo) = a— {tpN - V'oTv) 
Po 




(a - tto) 



(3.6) 



(3.7) 



y 



aBa dp 
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where a is is the half radius of the last closed flux surface (LCFS) at the height of 
the magnetic axis, Ba is used to normalise all magnetic quantities]^ = ip/a^Ba, 
po = p{i>o), qo = lii^o) and p is the GS2 flux surface label, which can have varying 
definitions, but which is always at the magnetic axis, increasing to 1 at the 
last closed flux surface (LCFS). In this work p is defined to be a^/a (the default 
choice), where a^{i') is the half radius of the flux surface at the height of the 
magnetic axis (note that a^{il)LCFs) = «)• 

The solution of the gyrokinetic equation depends on (among other quantities) 
the equilibrium gradients, dus/d-ip, dTg/dip and dco/d^. To simplify the treatment 
of these gradients, we expand the equilibrium quantities around the flux surface 
%jj = ipQ. Thus we may write: 



dng 
dip 

dT^ 

dip 

uj{ip) = u (ipo) + {^- i^o) ^ 



ris (ip) = Us (ipo) + - V'o) 
Ts{iP) = %{iPo) + {iP-iPo) 



^0 



V'o 



1 - 



sO 



X 

UJo + - 

a 



PqX 


dp 


a 


qo a 


dipN 


■ipNO _ 


PqX 


dp 


a 


qo a 


dipN 





dp 



dip 



(3.8) 

(3.9) 
(3.10) 



where we have defined the density and temperature gradient scale lengths, 

a 1 dn 



and 



ris dp 
1 dT, 



Lts Ts dp 
and the perpendicular flow shear 

qo dp 



PO 



PO 



(3.11) 



(3.12) 



po 



(3.13) 



where a^, the sign of ip, is defined in Section |33ln It will be immediately obvious 
to the reader that (owing to a historical quirk) 7^; is not in fact the perpendicu- 
lar flow shear, but is equal to the perpendicular flow shear multiplied by B/B^. 
However, in general B/ B^ ~ 1, so the parameter 7^; may reasonably be referred 
to as "the perpendicular flow shear". 



* Ba is defined to be the magnitude of the magnetic field at the centre of the LCFS. See the 
discussion of normalising quantities, Section l3.6.1l 

^ The introduction of this sign maintains the relation it • V = x^Ed/dy (see Section [3.4.21 
regardless of the sign of i/'/ which determines the sign of y. 
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3.2.2 Spectral Coordinates 



Since the gyrokinetic equation (in the absence of flow shear) has no explicit 
dependence on either x or y, it can be solved spectrally in both perpendicular 
directions. The reasons for doing this are firstly that the linearised gyrokinetic 
equation can be solved separately for each Fourier mode (see Section I4.1[) , sec- 
ondly that the pseudo-spectral calculation of the nonlinear term leads to faster 
convergence with perpendicular grid size (see Section I4.2|) and thirdly that the 
quadratic nonlinearity in the gyrokinetic equation is exactly energy conserving 
for a spectral system |77l,l78l] . A perturbed quantity A is therefore written as 



(3.14) 



where and ky are the actual perpendicular coordinates used by GS2 and we 
have defined the inverse Fourier transform ^ Historically, however, to make 
contact with the ballooning approximation, (|3.14|) was written ||6^: 



no, do 



where 



5* = no{a + q9o) . 
If we Taylor expand q, we can write 



S = no 



dq 

{a - ao) + - ipo) ^ 



+ no (ao + qoOo) 



(3.15) 



(3.16) 



(3.17) 



and if we absorb the term no (ao + qodo), which is constant over the simulation 
domain, into A, we may write 



^ ^ no dp 



Uq dp 



IpON 



a (a — ao) 



dip 



N 



a dipM 
kxX + kyU 



X 



dp 
riQ dp 



a dip 



N 



V'OiV 



(3.18) 



where ky = (no / a) dp/ dipN\^^^, = kysOo and s is the magnetic shear (see below). 
Thus, the expressions (|3.14|) and (|3.15|l are equivalent up to a constant factor]^ 



*My thanks to M. Barnes for his notes on this subject. 
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3.2.3 Magnetic Shear 

In the case of finite magnetic shear, that is, where § = {p/q)dq/dp > 0, the 
magnetic field lines are twisted between one flux surface and another, as is illus- 
trated by Fig. 13.21 which shows sets of magnetic field lines on two different flux 
surfaces (that is, at two different values of x). At finite magnetic shear the field 
line following coordinate system {tp,a,9) becomes non-orthogonal as q, and hence 
a, becomes a function of ip. Thus, the twisted shape of the magnetic field lines 
has been absorbed into the coordinate system. In the local coordinate system (x, 
y, 6), y becomes a function of x9, and as a result the full radial derivative of the 
eikonal dS/ dx becomes a function of 6. This has important consequences for the 
parallel boundary conditions, as is discussed in Section 14.1.41 




Figure 3.2: The shape of magnetic field lines on two neighbouring flux surfaces 
with (left) and without (right) magnetic shear. Each picture shows two planes 
of constant x within a flux tube, each of which encompasses a set of magnetic 
field lines. It can be seen that in the case of finite magnetic shear the gradient 
of any quantity taken between the two flux surfaces (that is, the effective radial 
gradient) can change dramatically along the field lines, even though the radial 
gradient in the local coordinates may be constant. 



3.3 A Note on Signs 

Although GS2 uses the {x, y, 9) coordinate system, at various points it is math- 
ematically convenient to use any one of 8 other coordinates: ip, a, 0, 9, r, p, R and 
Z, the vertical distance above the plane 9 = (the midplane). Although all of 
these coordinates have been defined above, the sign relation between all of them 
is not immediately clear: for example, in which directions, relative to the mag- 
netic field, do the toroidal and poloidal angles increase? 

We start by noting that there are two fundamental vectors to which all other 
vectors can be related: B^4^ , the toroidal magnetic field, and Ip, the plasma cur- 
rent which gives rise to the poloidal magnetic field We next define a set of signs 

^ In a stellarator, of course, there is no plasma current, and so aj does not exist. The choice of 
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Signs 



as = sign (/p ■ Bt) 
ai = sign {Ip ■ V0) 
(Tt = sign {B ■ V0) 
ae = sign (B ■ S/O) 
= sign {Vip ■ Vr) 
CTp = sign ( Vp ■ Vr) 



o-g = sign (g) 



= sign((V(/) x Vi?) ■ VZ) 



Table 3.2: Signs which relate the various coordinates used by GS2. 

(Table |3^ . The sign as is determined by the equilibrium field and is thus a given. 
Clearly aj = (Xto-B, and = ue by definition. Since also by definition 



where is determined by the magnetic field. In a system with up-down symme- 
try, az is, of course, irrelevant. Finally, although p can be defined in various ways 
(usually r /a) its always defined to be at the magnetic axis and 1 at the LCFS; 
thus, dp = 1, and = sign {yip ■ Vp). 

In a non-axisymmetric system, this is all that can be determined. In other 
words, given a magnetic field we can choose both at and ag, that is, we can choose 
the definitions of both (f) and 9. 

In an axisymmetric system, we can go further. Using Ampere's law, V x Bp = 
Ip , we can write 




(3.19) 



we find that aq = atae- For the sign of Z, we can write 



az = (Tt sign [{Bt x V-R) ■ VZ] = ata^, 



(3.20) 



sign [{B X /) ■ Vr] = 1, 



(3.21) 



which means that 

sign [{B X V0) ■ V^] = aja^. 
Using the usual expression B = Va x Vip, we can write 



(3.22) 



sign {\Vipf Va ■ V0) = cr/cr^, 



(3.23) 



which in an axisymmetric system means that 



aia^ = 1, 



(3.24) 



signs is much less restricted because the device is not axisymmetric, and there is no equivalent of 
the relation (I3l24t . 
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which means that 

a^ = (Tb, (3.25) 
ai = ag = = ataq. (3.26) 

In other words, we only have the freedom to choose the direction of either (p or 9. 

This restriction comes from the fact that in an axisymmetric up-down sym- 
metric system, all the equations must be invariant under a 180° rotation of the 
magnetic field about a horizontal axis which passes through the centre of the 
tokamak. Such a rotation would change aj and at: if we change the sign of (p to 
undo this, we must also change the sign of 9 to keep the equations the same. 

Lastly, we note that in Js^, where B = V-ip x Va (where ip = —ip), 

CTg = -ctb, (3.27) 
0"/ = -o"e = -o"^ = -^tCTq, (3.28) 



and so in Ref . |53l] , the sign of q is always opposite to this work, and the sign of 
either the toroidal field = RB ■ V0 or the poloidal field Bq = rB ■ must be 
opposite. 



3.4 Geometric Operators in the GS2 Coordinate 
System 



Certain of the operators listed in Table 13.11 can now be written explicitly in 
the (x, y, 9) coordinate system, specifically it ■ V, b ■ V (in certain circumstances), 

[h X V(<^)^J ■ Vhs, Vtfj-ibx V) and OFs/dtp. 

3.4.1 The Gradient Operator V 

In this section, we make heavy use of the gradient operator 



V = x|V.|| + !)|V,|| + Si,^ 



where x = Vx/ \Vx\, y = Vy/ \Vy\, 

\Vx\ = 

|V|/| = 

and 



a— \ViPn\ 
Po 



dip 



N 



dp 



iVal 



PO 



kp{9) = b-V9, 



(3.29) 

(3.30) 
(3.31) 



(3.32) 
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is not a wavenumber but is the parallel gradient (called gradpar within GS2) and 
is, in general, a function of 6*1^ 

3.4.2 The Equilibrium Flow Operator u ■ V 

In the case of an axisymmetric field, this operator, corresponding, to convec- 
tion by an equilibrium flow u = uR^Vcj), can be written very simplyO Using our 
general expression for the gradient operator 13.291 we can write 

x\Vx\ — +y\Vy\ — + hk,-^ (3.33) 

The third term in brackets is clearly 0, and the first term is as a result of axisym- 
metry. Using p.3|) and axisymmetry we can write 

V(t)-y = Vcj)- = T^^(t) ■ (V</' + qVa + OVq + Vz/) 

|Va| |Va| 

^ (3.34) 



|Va|i?2 



so that 



■u ■ V = 



dil>_ 



N 



dp 



PO 



d 

TT- (3-35) 

dy 



Using the local approximation (|3.10)) , then transforming to a frame rotating with 
velocity loq we can write 

d 

uV = a;7E — . (3.36) 

dy 

3.4.3 The Parallel Streaming Operator b ■ V 

In the GS2 coordinate system b ■ V can be written simply: 

b-V = fc,^. (3.37) 
The factor kp is calculated directly by the geometry module. 



® Except if the option equal_arc is specified, in which case gridpoints are evenly spaced along 
the field line and kp is thus a constant. 

' It is worth noting that it is only possible to have a strong (that is, |m| ^ vth) equilibrium flow 
where there is a high degree of axisymmetry 1,79.], and thus that the need for axisymmetry is not 
a restriction. 
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3.4. Geometric Operators in the GS2 Coordinate System 



3.4.4 The Nonlinear Term {c/B) (h x V((/?)^J • Vhs 

Using the expression for the gradient operator (|3.29[) , we may write 



- (6 X V(vp)^J ■ 



B sinv \ dy dx 



_ c \qo\ B 
B po Ba 



dip 



N 



dp 



dx dy 



PO 



\ dy dx dx dy 



(3.38) 
(3.39) 



where sinv = (x x y) ■ b (and thus |VV^| |Va| = Bsinv), and where we have 
ir te 

A/;, = ( -^^^ - -^^^ 1 , (3.40) 

(3.41) 



defined the nonlinear term, calculated in Section l42l as 

c (d(^dh_d(^dh 



where 



Ba \ dx dy 
dipN 



dy dx 



m 

Po 



dp 



PO 



and is called kxf ac within GS2. 

3.4.5 The Radial Perturbed Flow (c/B) {b x V{ip)^J V^j 

The radial component of the perturbed flow ■ (6 x V(v5)^ ) can now be 
written very simply. Using the general expression for the V operator, (|3.29[) , we 
can write 



-(bx V(¥^)^J- 



B sm V dy 



c 



dip 



dip 



N 



iVal d{ip) 



aBa 



dp 



dp 
d{^) 



po 



sin V dy 



po 



dy 



(3.42) 



3.4.6 The Background Gradients dFjOip 



Using the expression for the equilibrium distribution function (|2.7|) , we may 
write: 



dip 



d_ 

dip 

Fs 



risiip) 



27rTs{iP) 



3/2 



-es/Ts 



1 dus 1 
Us dip 




3^ 


1 dTs' 


a 


2y 


Ts dip_ 


1 dus 1 
ris dp 




3^ 


1 dTs' 




2y 


Ts dp\ 



dp 



(3.43) 
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Using the local approximation, (|3.11|) and (|3.12[) , we may write 



+ 



dp 



dip 



(3.44) 



3.5 The Gyrokinetic Equation in the GS2 Coordinate 
System 

Having written the geometric operators in the GS2 coordinate system we may 
now substitute them into the gyrokinetic equation (|2.10|) , giving 



f d d \ dh 



Rs 



+ 



T 

£s_ _ 3 
Ts 



2 L 



Ts 



Fsc d{ip) 



msliip) go 
TsB Vo J Baa dy 



Rs 



(3.45) 



where GS2 always makes the choice = 1. 

This is the equation that GS2 solves in an axisymmetric toroidal system. The 
limits of a cylinder and a slab can easily be derived from this equation. GS2 can 
also solve the gyrokinetic equation in a non-axisymmetric system, for example, a 
stellarator. This, however, is not covered in this work. 



3.6 Normalisations 

As is standard, normalisations within GS2 are chosen such that the nor- 
malised quantities are of order unity. First, a set of normalising quantities are 
chosen which are easily determined for a real-life system, and then all quantities 
within the gyrokinetic equation are normalised appropriately. Perturbed quan- 
tities are scaled up by a factor of egk to keep them, in the general case, of order 
unity. 



3.6.1 Normalising Quantities 

There are several things to bear in mind concerning the normalising quantities 

used in GS2, which are listed in Table l33l The first is that the reference species 

is not one of the species within the problem, for example ions or electrons, but is 

a hypothetical species to whose properties all the properties of the actual species 

are normalised. If, for example, the density or temperature of one of the actual 
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Quantity 


Definition 


a 


half diameter of the LCFS at the height of the 




magnetic axis 


T 


temperature of ref . species 




mass of ref. species 


rir 


density of ref. species 


Zr = l 


charge of ref. species 


Vthr 


^y2Tr/mr 


Ba 


toroidal magnetic field at the centre of the flux 




surface at the height of the magnetic axis 




ZreBa/rUrC 


Pr 





Table 3.3: Normalising Quantities: the definitions of the quantities used in the 
GS2 normalisations. 



species is chosen to be exactly 1.0, then the density or temperature of that species 
will happen to be the density or temperature of the reference species, but other- 
wise the density or temperature of the reference species will not correspond to 
any real quantity. However, the properties of the reference species can easily be 
calculated using the normalised properties of any actual species and the real val- 
ues of those properties for that actual species in the system of interest. Note also 
that the reference charge, Z^., is always equal to 1. Thus the normalised charge of 
electrons must always be set to -1, and that of hydrogenic ions to 1. 

The second is that the quantity a does not actually appear in the original 
gyrokinetic equation, or in its derivation (see Ref. |5^). Thus, although many 
quantities are normalised to it, and although it is given the definition "the half 
diameter of the LCFS at the height of the magnetic axis", in certain simplified 
geometries, most notably the slab, the cylinder, and shifted circle geometry, it 
has no physical significance. In more complicated geometries, its only physical 
significance is in the calculation of some of the geometric factors listed in Table 

m 

The third is that the normalising magnetic field Ba will in many simpler mag- 
netic geometries be simply equal to the magnetic field strength at the magnetic 
axis. Its rather involved definition becomes important in more complicated ge- 
ometries. 

3.6.2 Normalised Quantities 

The normalised quantities used within GS2 are laid out in Table l334l Table 
13.41 covers only the quantities and variables contained within the closed system 
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Name 


Quantity 
— y 


Definition 


GS2 Name 


Eqbm. Dist. Fn. 


Fns 


Fs/ins/vl,) 

/ \ ^/ zns/ 




Perturbed Dist. Fn. 


^Ns 


hs/{Pr/o)Fs 




Comp. Pert. Dist. Fn. 


9Ns 


qsl{prlo)Fs 


g 


Pert. Electric Potential 


7^ iv 


LP 1 ( Pr 1 Ci){Tr 1 e) 


phi 


Time 




tl ialvthr) 


time 


Parallel Velocity 


V\\N 


V\\IVih<, 

/ LI Lb 


vpa 


Perp. Velocity Squared 


N 




verp2 


Thermal Velocities 


VthNs 


Vfhtlvthr = \/Tj\[^/mi\[c 

LI Lb / LI LI W ly b / ly b 


stm 


Gradient Operator 


V/v 


aV 




Parallel Gradient 






gradpar 


Radial Coordinate 


Xn 


X/ Pr 




Poloidal Coordinate 


Vn 


V/Pr 




Radial Wavenumber 


kx N 




akx 


Poloidal Wavenumber 






akv 


Magnetic Field 


i\ 


BlBa 

j (1 


bmag 


Magnetic Flux 


I^N 






Current 


In 


I/aBa 




Nonlinear Term 


J^Ns 


UsIiclBaMTrleMFsIa^) 




Collision Operator 




CI {Vthrlci){Prlo)Fs 
1 \ LiLi 1 J v/^ / / y b 




Energy 


S AT a 

^ 1\ s 


6 a /To 

^ s / s 


enersv 


Magnetic Moment 


r^i\ s 


a/(TJ2B„) 




Lambda 


Xs 


PNs/^Ns 


1 


Flow Shear 


lEN 


lE/{Vthr/a) 


g_exb 


Temperature Gradients 


Kg 




tprim 


Density Gradients 


^ns 


O'l Lns 


f prim 


Temperatures 


Tns 




temp 


Densities 


nNs 


ris/rir 


dens 


Masses 


TUNs 


ms/rrir 


mass 


Charges 


Zns 


Zg/ Zr = Zg 


z 



Table 3.4: Normalised Quantities: definitions of the normalised quantities used 
within GS2, along with the names used for those quantities within GS2, where 
appropriate. 



of equations consisting of the gyrokinetic equation, the quasineutrality condition 
and Ampere's law. It does not cover the myriad diagnostic quantities which can 
be output from GS2, which are not covered in this work, since we are concerned 
with only two: Q and H, whose normalisations are given in Chapter ITU 

3.6.3 Normalisation of the Nonlinear Term 

Defining 
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Ns — ^2 



9{vn)r^ dhNs d{<^N)R^ dhNx 



dy^ dx 



N 



dx 



N 



dy 



N 



(3.46) 



we find that 



Ns 



(3.47) 



3.6.4 The Normalised Gyrokinetic Equation in the GS2 
Coordinate System 

Having defined a set of normalised quantities, we now substitute them into 
the gyrokinetic equation (|3.45|) to give 



Vthr Pr ^ 

a a \dtN 

Vthr Pr 



d d 

+ xnIen- 



dy 



N 



, , , Vthr Pr r.\ ^ths , dhg 
riNs + ( i's ] V\\NKpN 



a a 



+ 



\ C TrFs 

Fs I VdNs ■ Vtv^ATs + ^-A/tvs 

a a / Ba e 



^ Vthr 
Vthr Pr 



T, 



Ns 



Vthr Pr 



a a 



dt 



d d 

+ xnIen 



N 



dy 



a a 

(Vn) 



de 

Fs ) {Cn [hNs])R^ 



N 



Rs 



+ 



^Ns 2 ' 



msVthsVthr In qo \ c Tr F^ d{ipN) r^ 



Using the identity 



C Fa 



Bn Po 



1 Vthr pr 



Ba e a2 



dy 



(3.49) 



Ba e a? 2 a a ' 
and dividing by {vthr/a) {pr/ci)Fs, we obtain the normalised axisymmetric gyroki- 
netic equation in GS2 coordinates: 



d d \ dhj\i 1 

^7 \-XnJEN^ hNs+VthNsV\\NkpN—^ + VDNs-'^NhNs + --^Ns — {CN[hNs])R^ 

ot]\f (jyN J z 



d 

Ns \dtN 



+ 



l^ns \ ^Ns — I f^s 



+ xnIen 



0"^ 



d 



dy_ 



N 



{Vn)r^ 



2 In qo \ 1 d{ipN)R^ 

lEN \ ET-^- (3-50) 



VthNsBN Po '""J 2 dy 

3.6.5 The Normalised Quasineutrality Condition 

In a similar way, we can write the normalised quasineutrality condition: 

Z'^^^ipN = Zs I '^nNs{FNshNs)r- (3.51) 

"—^J-Ns „ J ^ths 
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3.7 Non-Trivial Geometric Operators 

The geometric operators in Table 13.11 which are not trivial to calculate in the 
GS2 coordinate system are calculated by the GS2 geometry module on a case- 
by-case basis, depending on the specific magnetic configuration chosen. These 
operators are kp, b x (b ■ V6), b x (Vlni?) and k±. In certain simplified geome- 
tries, for example, a slab, a cylinder or shifted circle geometry, they can be written 
as simple functions of 9, but in the general case, for example a numerical equi- 
librium taken from experiment, they must be calculated from the data for each 
grid point in the simulation. In this section, we describe how these operators are 
related to the normalised quantities that are calculated by the geometry module. 

3.7.1 The Parallel Gradient 

As noted earlier, the parallel streaming operator is written 

&-V = &-Ve^ = M^)^. (3.52) 

The parallel gradient kp, called gradpar within GS2, is calculated by the geometry 
module. 

3.7.2 The Magnetic Drifts 

The magnetic drifts are among those terms which are not rendered simple to 
express by the coordinate system used in GS2; hence, we do not expand Vds in 
the gyrokinetic equation (|3.50|| . 

For convenience of calculation, GS2 splits up the magnetic drifts into four 
components in the following way. First we use Equation (127) from Ref. |53^ : 

-jxB = Vp, (3.53) 

c 

where j = (c/47r)V x B and p is the total pressure, which expresses the balance 
between the magnetic field and the equilibrium pressure gradient (neglecting the 
high Mach terms), to write 

Bx{B- VB) = AttB xVp + Bx [{VB) ■ B] . (3.54) 

Using the identities B x {B ■ VB) = B^B x (b ■ Vb) and B x [{VB) ■ B] = 
B (B X VB) we can then write 
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1 



Vds = 7^ X vfib ■ Vb + -^VB 



2B 



53-11 



(3.55) 



Using (|3.15|) and normalising the magnetic drifts we may now write 



V, 



DNs ■ y N'l'Ns 



avthr 

Writing 



-Lr {v^iN + vIn) Bn X VnBn + ^%Bn X Vps 

AT -Dai 



N 



'N 



V NO — — r^yN^ — 
Pr op 



VAra + -VnQ , 



kyl\fS 



and using the identity 



we may write 

_ ih^skyN Ts 1 



IT. 1 



VnS. (3.56) 

(3.57) 
(3.58) 



Tr Zs 



2 di>_ 



N 



B% dp 



b X Vat-Bat 



VATtt + -, rVATg 



1 dij_ 



'N 



b X VN(3a ■ VTva + 



kyNS 



kyNS 

Vat? 



(3.59) 



"^5^ dp 

where Pa = 8np/Bl. Remembering that b x Vp ■ Vg = because the equilibrium 
pressure is a flux function, we may now write 

-Tns 



VdNs ■ Vat/IaTs 



where 



''±N 



2 [ 2 V 



kxN 



-UDNhNs 



kyNS 



+ '"WN \ ^kN + 



kyNS 



kN 



and where we have defined the four normalised magnetic drifts 

2 dipN 



(0) 



b X VnBn ■ VatQ; 





VBN 



{0) 



Bl dp 
2 djjN 

'^flNT^ ( b X VAT/3a ) ■ VatO + CO^BN, 



b X Vat-Bat 



5^ 



which are output by the GS2 geometry module. 



10 



Known as gbdrift, gbdriftO, cvdrift and cvdriftO respectively within GS2. 



(3.60) 
(3.61) 

(3.62) 

(3.63) 

(3.64) 
(3.65) 
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3.7.3 The Perpendicular Wavenumber 

The perpendicular wave number appears in the Bessel functions which arise 
as a result of gyroaveraging spectral quantities. Although the perpendicular 
wave number can be written simply in the GS2 coordinate system: 



kl = \VS\^ = kl \Wy\' + kl |Vx|^ + 2k^ky |Vx| \Wy\ cos 



V 



(3.66) 



the factors |Vx|, \Vy\ and cost; are in general functions of 6 and nontrivial to 
calculate. GS2 uses the definition (|3.16[) to write k\ in terms of three numbers: 



where 



kliO) = kl 



91 + 2-^fi'2 + 

KyS 



kyS 



93 



(3.67) 



91(0) 

92{e) 

93(0) 



|VAra| , 



(3.68) 
(3.69) 
(3.70) 



These three numbers are calculated by the GS 2 geometry module for the particu- 
lar geometry in use. It should be noted that g2 and are both equal to wherever 
s = 0. 



3.8 Geometric Factors in Specific Geometries 

3.8.1 A Note on the — a" Formalism 



As described by Miller |80l] . the "s — a" formalism |81|1 is a method for perturb- 
ing a general magnetic equilibrium by changing only two parameters, dp/d^', the 
gradient of the pressure, and dq/dip, the gradient of the magnetic safety factor. 
Changing the pressure gradient causes the magnetic flux surfaces to shift one 
way or the other according to the Grad-Shafranov equation; this is known as the 
Shafranov shift. Changing the safety factor gradient (i.e. the magnetic shear) 
changes the angle between the magnetic field lines in one flux surface and the 
next, which has a dramatic effect on the stability of certain modes (as will be 
noted later on in this work). The "s — a" formalism does not specify the magnetic 
equilibrium to be perturbed; hence, referring to "s — a" geometry is misleading. 

Two very commonly used geometric approximations which employ the "s — 



a" formalism are shifted circle geometry and Miller geometry |8C1|] . Shifted circle 
geometry, which is used in this work, is described in Section 13.8.21 
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3.8.2 Shifted Circles 

The shifted circle approximation describes a magnetic equilibrium composed 
of nested circular flux surfaces, generally, as here, taken in the limit of high aspect 
ratio. These flux surfaces may be shifted radially in the mid-plane by the Shafra- 
nov shift tts, and the field lines may be sheared from one surface to the next 
according to the magnetic shear parameter s. The geometric factors used by GS2 
can be written explicitly for shifted circle geometry and are given in Table 13. 5[ It 
should be noted that the exact geometric factors which are implied by a "shifted 
circle equilibrium", or an "s — a" equilibrium are not universally defined, and the 



different conventions can have significantly different results |82l ] (indeed, several 
different options are available within GS2). Here we merely quote the values 
used in this study. 



Quantity 


Value 


GS2 Name 




{1.0 + (r/R) cos9)~' 


bmag 


dp/dipN 


a/r 


drhodpsi 




1.0 


kxfac 


^kN 


(2a /R) (cos 9 + {s9-as sin 9) sin 9) 


cvdrif t 


/ 

^kN 


{2a/R)ssm9 


cvdrif to 


^SJBN 


{2a/ R) (cos 9 + {s9- a, sin 9) sin 9) 


gbdrif t 




{2a/R)ssm9 


gbdriftO 


9\ 


1 + (s sin 6* — as sin 9)'^ 


gds2 


92 


— s {s9 — as sin 9) 


gds21 


93 




gds22 


kp 


a/Rq 


gradpar 



Table 3.5: Geometric Factors in Shifted Circle Geometry 



3.8.3 The Slab 

The slab is perhaps the most simple magnetic geometry possible, consisting of 
straight magnetic field lines with a constant shear. It can be derived as a reduction 
of the shifted circle equilibrium, by setting = 0, and letting r — )■ oo, i? — > oo 
but r/R — )• 0. The geometric factors used by GS2 can be written very simply for 
slab geometry and are given in Table 13.61 



3.8.4 The Cyclone Base Case 

The Cyclone Base Case is a special case of shifted circle geometry which has 
been used for benchmarks between gyrokinetic codes (see e.g. IZSt]). In effect, it 
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3.8. Geometric Factors in Specific Geometries 



Quantity 


Value 


GS2 Name 


— i V 


1.0 


bmae; 

o 


dpIdiJjN 

r 1 T i V 


1.0 


drhodpsi 




1.0 


kxf ac 


^kN 





cvdrif t 


1 





cvdrif to 







gbdrif t 







gbdriftO 


91 


1 + (s sin 6 




92 


-s^e 


gds21 


93 




gds22 


kp 


kp 


gradpar 



Table 3.6: Geometric Factors in Slab Geometry. The quantity kp is an input pa- 
rameter to the simulation. 

is shifted circle geometry with the following choice of parameters: 



r 
R 


= 0.18, 


(3.71) 


a 


= R, 


(3.72) 


Q 


= 1.4, 


(3.73) 


Ki 


= 6.9, 


(3.74) 




= 2.2, 


(3.75) 


as 


= 0, 


(3.76) 


s 


= 0.8. 


(3.77) 



All simulations in this work use the Cyclone Base Case (with s = 0) as a starting 
point, with the exception of Chapter [131 which considers the case of the slab. 
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Chapter 4 



The GS2 Algorithm 



4.1 The Linear Implicit Algorithm 

4.1.1 Equations 

The linear solver advances the closed system consisting of the linearised 



gyrokinetic equation (|3.50|) for each species and the quasineutrality condition 
(133ID0 Starting from (|3.50|) and (|3.51|) and dropping the nonlinear and perpen- 
dicular flow shear terms, the equations to be solved are 



dh 



Ns 



dtN 



+ VthNsV\\NkpN 



dhs 



+ Vz 



Tns 

and 



dt 



N 



DNs ■ 



Vat/iats — (Cat [hNs]) 



2 In go 



Ns 



O -lEN I „ 
VthNs -DiV PO I ^ 



d?v /TP , 



ths 



1 d{^N) 



Rs 



dy 



We start by introducing the complimentary distribution function]^ 



9ns 



h 



Ns 



T, 



Ns 



■{^Ns) 



Rs' 



(4.1) 



(4.2) 



(4.3) 



^ As well as the parallel and perpendicular components of Ampere's Law in the electromag- 
netic case. 

^ This is done for historical reasons, although the complimentary distribution has the pleasing 
property that g — Ois the exact solution of the long-wavelength shear- Alfven- wave equation [83]. 



4.1. The Linear Implicit Algorithm 



for which the evolution equations are 



i 



— 1- VthNsV\\NkpN-^ \- VdNs ■ '^NQNs — {Cn [hNs])ii^ 



dt 



N 



T, 



VthNsV\\NkpN + VdNs ■ Viv(v?)^^ 



+ 



Ns 



de 



2 In qo \ l^(<^Jv)fi, 

-lEN 



VthNs Bn Po 



2 dy 



(4.4) 



and 



Ns 



U 



ths 



T, 



Ns 



I r 



(4.5) 



We now Fourier expand all perturbed quanties in the perpendicular plane 



using (|3.14)) to give: 



-7^ H VthNsVml^pN^^ H i-^(^DNgNs — {'^ N [iINs\) 

01 N 00 Zg 



T, 



Ns 



dJo (as) Tns t f 
VthNsV\\NkpN ^ h 1^-ujdnJo [as) V 



+ ^■ 



f^ns \ ^Ns 2 ' 



2 /jv go 

VthNs Bn Po 



lEN \ Jo (as) <^Af, (4.6) 



and 



Y] Zg'^ipN (1 - r,) = V [ ^nNsFNs{gNs)r- (4-7) 
-^Ns ^ J Vths 

where Jq is a Bessel function of the first kind and results from gyroaveraging a 
Fourier-transformed quantity, = y/T^s'Tn- NsV±Nk±N /ZsB^ = v±k±/D.s, and the 
integrated Bessel function 



^si9) 



—FnsJq (as) = exp 



-'ths 



'^±P. 



(4.8) 



where Jq is a modified Bessel function of the first kincj^. 

We proceed by collecting derivatives of each fluctuating quantity in (|4.4)) , 
dropping circumflexes and writing: 



dg 



Ns 



dt 



+ 



dgNs , , dJo (as) (fN , , dJo (a^) fN,,j, . 

+ O-OgNs = bt h bg — h 60 Jo [as) ^N 



de 



dt 



de 



(4.9) 



^Note that the collision operator must still be calculated in terms of hs 
* The quantity tinsZITs/Tns is called gamtot in the code. 
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where 

ae 

Co 

bo 



VthNsV\\NkpN) 
Tns 

i'—p;-^DN, 



0, 



T. 



Ns 



"VthNs^WNkpNi 



-iuDN + i- 



k 



yN 



2 In go 

VthNs Bn Po 



lEN 



(4.10) 
(4.11) 
(4.12) 
(4.13) 

(4.14) 



For the rest of Chapter H] we drop the subscript and work only in normalised 
quantities. 

4.1.2 Discretisation and Parallelisation 

The distribution function gs{9,kx,ky,e, \,a) is stored on a six-dimensional 
grid of size ntheta x nakx x naky x x nA x 2. The perturbed fields are functions 
only of 6*, k^ and ky. The linearised gyrokinetic equation (|4.9|) is in fact a ID partial 
differential equation in 9, with the values of the other coordinates appearing only 
as parameters; this means that its solution may be calculated separately for each 
value of kj;, ky, e and A, enabling very efficient parallelisation. The field equation 
(|4.7[) requires integrating over velocity and species, but can still be solved sep- 
arately for each value of k^ and ky. Thus, a typical strategy when solving the 
linear problem would be to divide it among nakx x naky x ne * nA * 2 processes: 
each would solve (|4.9|) independently, with communication only required for the 
solution of (l47b . 

The gyrokinetic equation (|4.9[) is discretized using a compact stencil to keep 
its inversion matrix bidiagonal. The calculation of the differences is complicated 
by the introduction of space centering and temporal implicitness parameters, tq 
and Tt, known as bakdif and f exp in the code. The space centering parameter 
re ranges from (centered) to 1 (a first order upwind scheme), and the temporal 
implicitness parameter varies from (fully implicit) to 1 (fully explicit). With 
these parameters, the perturbed potential and its differences may be written 



dip 


1 


'dt 


2 


dip 




dz 


= rt 




1 




2 



:i - re. 



n+l 



At 



+ [l + re) 



At 



+ (1 - n) 



n+l 



AO ' AO ' 

+ (1 + re) [r,^r+i + (1 - n) } 



(4.15) 

(4.16) 
(4.17) 
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and similarly for the distribution function. The special choice of = and 
rt = 1/2 corresponds to a space-centered, time-centered (trapezoidal) scheme, 
first proposed by Beam and Warming l Is^ P which is second order accurate and 
unconditionally stable. 

Substituting the differences into (|4.9|) , we obtain the discretized gyrokinetic 
equation (suppressing the species index s and defining (p = Jq^)'- 



n+l _ n „n+l _ n 



At 



+ ae 



At 

rt TT^ ^ (1 - rt 



A9 



A9 



+ f { (1 - [rt9? + (1 - gr'] +{l + re) [rt9?^^ + {1 - n) ^I^Y] } 



rt TT^ ^ (1 - rt) - 



A9 



A9 



+ j{a-ro) + (1 - rt) ^r'] + (1 + re) [rt^7+i + (1 - n) ^Iti] } ^ (4-18) 

which may be written as 



A,g^ + A,g^^, + B,g^+' + B,g^^,' = D,^^ + D^^^+i + E.^T' + ^2^r+Y (4-19) 



~n+l 



~n+l 



where 



B2 



I - re rt ao 



2At 

I - re 

2At 

l + rg 



^ — rt ao 
- ^s^Q- + y (1 - re) (1 - n) , 



(4.20) 
(4.21) 
(4.22) 
(4.23) 



The quantities 5^^, B2^Bi, A2 and Ai are called ainv, r, a and b in the dist_f n 
module of GS2. The definitions of Di, D2, Ei and E2 are unimportant as the RHS 
of (|4.19|) is calculated in a different way within the code (see Section |4.1.5|| . 



^Not to be confused with the Beam-Warming scheme, a different scheme from the same au- 
thors. 
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4.1.3 Advancing Implicitly using Green's Functions 



In each iteration the new distribution function must be calculated using (|4.19|) 
(reproduced here for convenience), and the new potential via the discretized 
quasineutrality condition 14.241 



E Z'^Usjl Ts) ^ +1 _ , , n+1 

s ^ s,£,X,a 

where JT^.a.ct are the weights for the velocity integral and where we have defined 
the velocity and species integration operator V. Since this is an implicit scheme, 
solution of (|4.19|) requires <^"^^ and solution of (|4.24[) requires g'^^^. Generically, 



this can be cast as a single matrix equation, writing (|4.19|) as 



and (I4.24|l as 



A ■ g^^^ + R-9'' = R- + 1 ■ <^", (4.24) 

F ■ (^"+^ = G ■ ^"+1 (4.25) 



so that 



^"+1 = • ^" - ^ • (4.26) 

where A — F_ are square matrices of size Ntot = ntheta x nspec x ne: * nA * 2 and 

K= {A- RJoEr^0 (4.27) 

is the implicit response matrix. However, initial inversion of the (dense) response 
matrix would require O (Nf^^) operations, and subsequent timesteps would re- 
quire O (A^toi) operations: for any reasonable problem size the computational cost 
quickly becomes prohibitive. 



Kotschenreuther et. al. |67l] realised that this cost could be avoided as follows. 
Since the distribution function responds linearly to changes in the perturbed po- 
tential, the response of the distribution function to any arbitrary change in the 
perturbed potential can be calcuated as a superposition of the responses to a com- 
plete set of delta function changes to the perturbed potential: a Green's function 
approach. This allows the new field to be calculated before the new distribution 



function, after which the new distribution function can be found using (|4.19|) 
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The procedure is this. First define ip* as the change in the potential from one 
timestep to the next: 

^* = _ (4 28) 

Next split g"^^^ into two parts: 



(4.29) 



and define g'' as that part of the new distribution function which results from the 
change in the potential ip*: 



B.g'r + B.gt^.^E^ + E^^,. 



(4.30) 



Substituting (14281) . (I429b and (l430b into (I4l9b we find that is that part of the 



new distribution function which results from the old distribution function and 
the old potential]^ 

B,g^ + B^gU = -Mgl - A^g^^, + D,^"^ + D^^^+i + ^i^^ + ^2^IVi- (431) 



The motivation for the definitions (|4.28|) , (|4.29|) and (|4.30|) lies in the quasineutral- 
ity condition. Firstly let us define the response matrix, the response of the new 
distribution function at i to a change in the potential at j, 



^g \ 

^ J = {Bi6i^k + -82^2,^+1) {EiSk,j + E25kj+i) 



by rewriting (|4.30|) as 



(4.32) 



(4.33) 



Secondly let us substitute the definitions (|4.28|) , (|4.29|) and (|4.33|) into the 
quasineutrality condition (|4.24[) to give 



which can be rearranged to give 



p 



(4.34) 



(4.35) 



^ The superscripts p and c, which stand for predictor and corrector, are used here because the 
algorithm resembles a predictor-corrector scheme. However, as show by Belli [69], the linear col- 
lisionless algorithm presented here is in fact mathematically equivalent to the implicit algorithm 
presented in Ref . 16711 
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where the matrix is a square matrix of size ntheta, which must be calcu- 
lated once at the beginning of a simulation. Calculating the matrix is really a 
matter of calculating the response matrix {5g/5(p)^-. This could be done by direct 
calculation from i?2, Ei and E2, but this effort can be spared by noting that if 
V?* is set equal to 5jk then g^, the solution of (|4.33|) , i.e., the solution of (|4.30|) , is 



actually the /cth column of the response matrix. 

The algorithm for GS2 may now be written very simply indeed: 

1. Calculate each column of the response matrix by setting ip* = 5jk and solv- 
ing (I4.30D for k= 1, 2. ..ntheta. 



2. Calculate the matrix 

3. Specify the initial values of g and (p. 

4. Solve the equation (|4.31|) to find g^. 

5. Solve the equation (|4.35|) to find p* and hence (/j""*"^. 

6. Solve the full gyrokinetic equation (|4.19|) to get (7""*"^. 



7. Repeat steps 2-4 for every timestep. 

The effort of implementing this algorithm can be further reduced by noting that 
solving (|4.30|) is the same as solving (|4.19|) with g"' and (p"' set to zero, and (/?"+^ set 



to ip*, and solving (|4.31|) is the same as solving (|4.19|) with y}""*"^ set equal to p"'. 
And so now the algorithm may be written in a form that directly corresponds to 



the actual statements and subroutine calls in the code 



i 



With apologies for the corruption of the mathematical symbols by the use of the imperative 
programming style. 
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1. Set k = l. 

2. Set if* = 5jk, ¥5"+^ = ^" = and = 0. 

3. Solve the equation (I4l9ll for 51"+^ 

4. Set (5^7/5^),, = ^ 

5. Repeat steps |2]-S1 for k = 2...ntheta. 

6. Calculate the matrix P ^ 

7. Specify the initial values of (7" and (y?"^. 

8. Set 93"+! = 

9. Solve the equation (|4.19|) for 
10. Set = c/"+i. 



11. Solve the equation (|4.35|) to find 

12. Set = + 

13. Solve the equation (14191) to get ^"+^. 

14. Repeat steps [8l-[T3l for every timestep. 

Thus, there are only three major steps that need to be implemented: solving the 
gyrokinetic equation (|4.19|) , solving the equation (|4.35|) and calculating P^. The 



latter two operations require the application of the velocity and species integra- 
tion operator V (see Chapter [S]) and some linear algebra. Solving the gyroki- 
netic equation itself requires careful treatment of the boundary conditions, as 
described in the next section. 

4.1.4 Solving the Discretized Gyrokinetic Equation 



Solving (|4.19|) requires the correct treatment of the parallel boundary condi- 



tions. The cases of trapped and passing particles are treated separately. The 
boundary conditions may be summarized as follows: 

• For passing particles when the magnetic shear is zero, the boundary condi- 
tions are periodic (what leaves at one end must come back in at the other). 

• For passing particles when the magnetic shear is finite, the boundary con- 
ditions are zero incoming. 
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• For trapped particles, the distribution function for particles arriving at a 
bounce point must equal the distribution function for particles leaving it. 

The case of passing particles is complicated by the fact that in the case of finite 
magnetic shear s, the total radial derivative of the eikonal dS/dx becomes a func- 
tion of ^ H: 



Kt = K- kys9 (4.36) 



dS 
dx 

where kx is the "actual" radial wavenumber of a mode as defined in Section l3.2.2l^ 
The practical upshot of this is that while when magnetic shear is equal to zero the 
boundary condition is a very simple periodical one: 

g {kx, ky, 9 = 9max) — 9 {kx, ky, 9 = 9max) , (4.37) 

the boundary condition when the magnetic shear is finite is such that the end of 
one field line with a given kx connects to the beginning of another field line with 
the correct value of kx such that the value ofkxt matches at the join: 

g {kx, ky, 9 = 9-uiax) — 9 {^x ~ '^kyS9,„iax, ^y, 9 = —9max) , (4.38) 

where '2.9max is the length of the parallel domain (equal to vr x nperiod in the code). 
This is known as the twist and shift or linked boundary condition. Implementing 
it requires that the quantity 2kys9max/^kx (where Akx is the spacing in kx on the 
grid) be an integer: in the code that integer is called j twist. 

Thus, the set of nakx radial modes are linked into a set of nakx/ j twist parallel 
domains (see Fig. l4.1[) . The gyrokinetic equation (|4.19|) is solved for each parallel 
domain separately. Linearly, each connected parallel domain corresponds to a 
single ballooning modeE the boundary conditions for the connected domain are 
thus the same as those for a single ballooning mode: the distribution function 



must vanish at either end |67|1 . Since (|4.19|) is only first order in 9, this can only be 
enforced at one point in theta: physically this done by setting 9 = for incoming 
particles at either end of the domain: for particles with a > 0, g = for the lowest 
value of the extended theta grid, for particles with a < 0, g = for the highest 
value of the extended theta grid. 



It should be noted at this point that Ref |66|] defines as kys9a, which means that confus- 
ingly in GS2 the "actual" values of the radial wavenumber are sometimes referred to as values of 

' Note that if one wanted to simulate a single mode linearly, one could abandon this twist and 
shift arrangement and extend the parallel domain sufficiently far for convergence by increasing 
nperiod. This is the situation described in Ref. [67]. However, the twist and shift system is nec- 
essary for nonlinear runs, as it is necessary to have a large number of k^s to resolve the turbulent 
spectrum. 
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Figure 4.1: An example of the extended theta grid where j twist = 2. The crosses 
show the gridpoints and the lines show the connections between them. The nakx 
separate segments are connected into j twist domains. The left panel shows the 
domains vs. 9 and k^, which index the gridpoints. A parallel segment with = 
nAkx connects at = Omax to a segment with k^ = {n — j twist) Ak^ (see (|4.38|) ). 
The right panel shows the domains vs. 9 and k^t- The quantity k^t varies with 9 
(see (|4.36|) ). At either end a parallel segment connects to the parallel segment that 
has the same k^t at the opposite end. 



For trapped particles, the situation is much simpler, and can be written as 

g{a = +l,9 = 9^i,) = gi(T = -l,9 = 9ub), 

g{a = +l,9 = 9i,) = g{a = -l,9 = 9i,), (4.39) 

where_^„fe, 9ib are the points where the trapped particles are reflected by the mirror 
effect. 



10 



To recapitulate, the boundary conditions must be solved differently for three 
different cases. For the periodic case, the distribution function at one end must 
match the distribution at the other. For the twist and shift case, the distribution 



^"Since the magnetic moment /i = mv\/B of any particle is conserved, as a particle moves 
from the outboard side of a tokamak to the inboard, the increase in B causes an increase in v±. 
Since total energy e = m + of the particle is also conserved, its parallel velocity decreases 
until it reaches zero and changes sign — the particle boimces back. 
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function at one end of one parallel segment must match the other end of the con- 
nected parallel segment. At each end of the chain of connected parallel segments, 
zero incoming boundary conditions are applied. For trapped particles, the distri- 
bution of particles arriving at each bounce point must match that of those leaving 
it. The exact routine for implementing these conditions for the twist and shift case 
is too tortuous to be covered here, but we will describe the general principle and 
the specific cases of trapped particles and passing particles in a periodic system. 

The general principle follows the standard method for matching boundary 
conditions for a partial differential equation: first divide the equation (|4.19|) into 
a homogeneous and an inhomogeneous part, solve each part independently, and 
then combine them in a way that satisfies the boundary conditions. Both parts of 
the equation are solved separately for each sign of velocity, so that there are four 
solutions for each parallel segment. For the periodic case and the case of trapped 
particles, these four solutions can then be combined to satisfy the boundary con- 
ditions of each segment separately. For the twist and shift case, satisfying the 
boundary conditions requires communication between each segment of the con- 
nected parallel domain, which leads to significant complications (as exemplified 
by the routine init_connectedJoc). 

The equation for the homogeneous part g°'^'^ is 

B,gr + B2g'^ll = 0, (4.40) 
and the equation for the inhomogeneous part (7"^"" is 
B,gr + B,g^;r = -A,g? - A,g^^, + D,^- + D.^l,, + E.^^' + E.^^^l (4.41) 



All terms on the RHS of (|4.41[) are grouped together and called the source. Both 
(|4.40|) and (|4.41|) are solved first for a = -1 (f y < 0) and then for a = +1 {v\\ > 0). 

For passing particles, when a = —1, the initial conditions are g^tgiid = 1' 
fi'ntgiid = 0/ where the 9 index runs from — ntgrid to ntgrid and 

ntgrid = (2 x nperiod — 1) x ntheta. (4.42) 

Both equations are then solved starting at the highest theta and sweeping to the 

left. When a = +1, the initial conditions are g""^tgrid = h ^-ntgrid = 0- Both 

equations are then solved starting at the lowest theta and sweeping to the right. 

This results in the generic solutions illustrated in Fig. I4.2l fa). 

For trapped particles, when a = —1, the initial conditions are g^T-i,uh = 1 ^^id 

g'^T-iv.h = 0, where ub is the index of the upper bounce point. Both equations 

are then solved starting at O^b and sweeping to the left. When cr = +1, the initial 
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conditions are g°^^ 



+l,lb 



-1,1b 



-1,1b 



, where lb is the index 



of the lower bounce point. This ensures that the (|4.39|) is satisfied at the lower 
bounce point for any linear combination of g""''^ and g''-^'" . Both equations are 
then solved starting at and sweeping to the right. This results in the generic 
solutions illustrated in Fig. l4.2r b). 




max 







e 



ub 



9 



max 



Figure 4.2: Generic solutions of (|4.40|) (red, solid) and (|4.41[) (green, dashed) for 
(a) passing particles and (b) trapped particles. The arrows indicate the sign of v.. . 



For trapped particles and for passing particles in a periodic system, the four 



solutions may now be combined to satisfy the boundary conditions (|4.39|) and 
(|4.37|) respectively. Defining the total solution as 



5";ii = &+/5+&, 

g:t\=9T=., + P-9:T-i, (4.43) 

then for trapped particles the boundary condition at the lower bounce point re- 
quires /3+ = and the boundary condition at the upper bounce point requires 
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that: 

9a^-l,ub + f^-9a^~l,ub = 5'ct=+1,-u6 + f^+9a^+l,uby (4.44) 

which means that 

,new 

9a=+l,ub 

remembering that ^'""'Li^ub = 1 ^rid 9'^T-i^nh = 0- For passing particles in a peri- 
odic system, we require: 

ne-u) I o one _ new , o one 

y(T=±l, ntgrid K±yo-=±l,ntgrid ~ yo-=±l, -ntgrid /^±yo-=±l, -ntgrid' 

,new ^new 

,^ n "(T=itl, ntgrid "o-=±l, — ntgrid 46) 

„one „one ' \ • ) 

ytT=itl,— ntgrid "^=±1, ntgrid 

In a similar manner, in the twist and shift case, a set of linear simultaneous equa- 
tions can be constructed for the values of /3+ and /3 from every segment of a con- 
nected domain to satisfy (|4.38|) at each connection and the zero incoming bound- 
ary at either end. 

4.1.5 Calculating the Source Term 

The only complication in calculating the source term (that is, the RHS of (|4.41[) ) 



is the fact that the subcripts i and i + 1 switch with a — depending on which way 
along the field line the equation is being solved. To reduce coding effort, the RHS 
of (|4.41[) (which is actually the RHS of (|4.18|) plus the terms involving g"^), is tidied 



up in the following way. Defining 

(^- = ri<^" + (l-r,)<^"+^ (4.47) 

and 

(p^ = - (4.48) 



then the RHS of (144111 is equal to 



-A,g^ - A^gl^, + he^ + 'i [(1 - r,) + (1 + re) vZr] cj = +1 
-A,g^ - A,gf_,, + + | [(1 + r,) + (1 - re) ^^^i] ^ = -1- 



(4.49) 



Defining 

^' = r^e^tl^ + rT^r, (4-50) 
where the parameters and ^ are defined in the following way: 



11 



Called bdf ac_p and bdf ac_m in the code. 
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rl = l + are, 



r™ = 1 - are, (4.51) 
the expression (|4.49|) can be written as 



Ae 



(4.52) 



4.2 The Addition of the Nonlinear Terms 

4,2.1 Time Differencing 

In contrast to the linear parallel convection term v\\b ■ Vhs, the nonlinear term 
Af = -^b X V((y9)^ ■ Vhs is treated explicitly, and a Courant-Friedrichs-Lewy (cfl) 
condition is used to keep the time step small enough for stability. Accuracy is im- 
proved by the use of a linear multistep method: the third order Adams-Bashf orth 
formula. Thus, the nonlinear contribution to the discretized gyrokinetic equation 
(|4.19|) , which appears purely on the RHS and is added to the source calculated in 



Section 14.1.51 can be written [72(1 



|aC - Ik-' + ^K-'- (4.53) 

Where Afc is the discrete Fourier transform of A/" in the perpendicular directions. 
It remains to calculate Af. 



4.2.2 Calculating the Nonlinear Term 

In the field line-following coordinates defined in Section 13.21 the nonlinear 
term may be written (see (|3.38)) ): 



Af=nJ^^-^^]. (4.54) 
^ \ dx dy dy dx J 

The nonlinear term could be expressed and calculated as a complicated sum of 
the Fourier components of and g. However, it requires fewer operations IZZD to 
transform the fields and the distribution function back into real space, calculate 
J\f using (|4.54[) and then Fourier transform it to determine A^. 



Thus, the four partial derivatives in (|4.54[) can be written: 
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dh 

dy 



dh 
dx 



dx 



kyg + — - — ip 



T 



dy 



k^g H — ip 



(4.55) 
(4.56) 

(4.57) 

(4.58) 
(4.59) 



temporarily restoring circumflexes to denote spectral quantities and denoting the 
inverse Fourier transform defined in (|3.14|) as J-"^. This pseudo-spectral method 
of calculating the nonlinear term leads to spectral accuracy for the evaluation of 
the derivatives IZZD- 

There is an important complication to this procedure which may be sum- 
marised as follows: energy from the truncated set of Fourier modes used in a 
given simulation would, had a larger number of modes been used, have been 
scattered via nonlinear interaction into modes with wavenumbers up to 1.5 times 
as large as the highest wavenumber in the truncated set. If the above calcula- 
tion were to be done using only the truncated set, that energy will be artificially 
added to the truncated set as a result of aliasing. This effect can be removed by 
performing both the inverse and forward Fourier transforms in the calculation 
using a larger grid whose highest wavenumber (in both the kx and ky directions) 
is greater than 1.5 times the highest wavenumber on the actual grid. Before the 
inverse transform to get the real space fields and distribution function the extra 
gridpoints are set to zero; after the forward transform of the calculated nonlinear 
term the energy in the extra wavenumbers is simply discarded. 
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Chapter 5 

Velocity Space and The Collision Operator 



This chapter is a summary of Chapter 4 ofRef. h83il and ofRefs. h7Ci\7l\] . 



5.1 Introduction 

The distribution function is a field in a five-dimensional parameter space, and 
the gyrokinetic equation must be solved in five dimensions. It is of paramount 
importance to use the minimum possible resolution in every dimension (while 
still correctly resolving the physics) in order to keep the computational and mem- 
ory requirements of the simulation within the limits of what is attainable. The 
minimum resolution in real space is set by several factors, including the need to 
resolve both the driving scale and the dissipation scale. The minimum resolu- 
tion in velocity space is set by the need to resolve the small-scale structures that 
develop in velocity space as a result of processes such as Landau damping and 



phase mixing |86l,l83|1 



This minimum can be reduced in two ways. Firstly, we repeat the observa- 
tion that the gyrokinetic equation is solved separately for each velocity, and that 
the details of velocity space are important only for the velocity integration oper- 
ator V, and the collision operator C. The resolution requirements for the integra- 
tion operator can therefore be reduced by using spectrally accurate integration 
methods (although this complicates the implementation of the collision operator 
slightly). Secondly, we note that most systems have a finite coUisionality, which 
will lead to dissipation. The introduction of a collision operator which models 
the operation of these collisions will cause a smoothing in velocity space, and a 
reduction of the minimum velocity resolution required. Since the rate of dissipa- 



tion of turbulent energy is controlled only by the rate of its injection |86|l . and not 



5.2. The Velocity Space Grid 



by the coUisionality, it is possible to have a coUisionality slightly larger than that 
observed physically, leading to reduced velocity space structure and a further re- 
duced minimum resolution, without affecting results such as the calculated heat 
flux. 

As well as smoothing velocity space structure, collisions are also responsible 
for the irreversible heating of the background equilibrium that results from the 
turbulence ||86l. Since this heating must be included accurately in any attempt 
to attain a statistically steady state, it is desirable for the collisions to conserve en- 
ergy, particles and momentum and to satisfy Boltzmann's H-Theorum [[toI] . GS2 
now contains a collision operator which satisfies all these criteria arid provides 
diffusion in both parallel and perpendicular velocity dimensions [[toI,!/]! . 

This chapter describes the layout of the velocity space grid and the imple- 
mentation of the collision operator as they currently stand, without justifying the 
choices that were made therein. 



5.2 The Velocity Space Grid 

Since the perturbed distribution function is independent of the gyrophase, 
only two velocity space coordinates are required. In GS2 we choose the kinetic 
energ}|l] e and \ = ^/e (dropping the species subscript s), a choice which elimi- 
nates all velocity derivatives from the coUisionless gyrokinetic equation |85(1 . The 
layout of the numerical grids of these two coordinates is chosen to result in a 
spectrally accurate velocity integration operator. 

The energy integration operator can be written 

Vth / dx x^, (5.1) 
Jo 

where x = v/vth- The energy integral is split into two sections, from to xq and 
from Xq to infinity, where xq is a free parameter, usually chosen to be > 2.5. In 
the first section, the grid locations are chosen using Gauss-Legendre quadrature 
rules. The integral in the second section is performed by a change of variable to 
y = x^ — Xq, so that the integral becomes 

1 



dy e 



eyJy + xlG{x) 



(5.2) 



where G {x)'\s the function to be integrated. In this interval, the grid locations are 

chosen using Gauss-Laguerre quadrature rules. 

^It should always be remembered that the choice of energy as a velocity coordinate necessi- 
tates the addition of another index to the distribution function: the sign of the parallel velocity, 
a. 
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The choice of the grid layout for A is complicated by the existence of parti- 
cles trapped within magnetic wells in certain curved magnetic geometries. We 
require there to be a value of A for each parallel grid point, corresponding to the 
particle which bounces at that particular grid point. We also require there to be 
a concentration of grid points around the trapped /passing boundary in A — e 
space. 

Accordingly, the A grid is divided into two regions: one for passing particles 
with < A < 1/ Bmax and one for trapped particles with 1 / B^ax < A < 1/ Bmin- 
In the first region, the integration variable ^ = Vl — XBmax is chosen, and Gauss- 
Legendre quadrature rules are used to obtain the locations of the grid points. In 
the second region, values of A are chosen to correspond to the bounce points: for 
every value of 9, a value of A is chosen such that f y (9, A) = o|^ 

An important observation must be made here, namely, that since there must 
be a value of A for every bounce point, the user has only limited freedom in 
choosing the total number of values of A. Specifically, the user may only specify 
the number of untrapped values of A; the input parameter ngauss is equal to 
half this number. For example, with ntheta = 14 and concentric circular flux 
surfaces, there are 6 bounce points corresponding to different values of A, so nA = 
6 + 2 X ngauss. The number of bounce points depends on geometry, but is usually 
around half ntheta. 

The velocity grid for the case of ntheta = 14, ngauss = 4, negrid = ne = 8 and 
Cyclone Base Case geometry (concentric circular flux surfaces with r/R = 0.18 
and q = 1.4, see Section |3. 8. 4|) is shown in Fig. 15.11 



5.3 The Collision Operator 
5.3.1 Theory 



The collision operator defined in Ref . |70tl causes diffusion in both the energy 



and A directions, while locally conserving particles, energy and momentum. In 



its gyroaveraged spectral form, used in GS2 |7l|l , it can be written as the sum of 
five parts: 

CcKlh] = L[h] + D[h] + UL[h] + Uoih] + E[h], (5.3) 

where 



^ Note that in a stellarator with multiple magnetic wells, a single value of A may in fact corre- 
spond to several bounce points. 
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Figure 5.1: The layout of the velocity space grid corresponding to ntheta = 14, 
ngauss = 4 (corresponding to 14 values of A), and negrid = 8 (i.e. 8 values of e) in 
(a) E and A space and (b) v± and \v\i | space at the outboard midplane (B = Bmin)- 
It should be remembered that there are twice as many velocity grid points as 
displayed owing to the two signs of v\\, indexed by a. 



(5.4) 

is the gyroaveraged Lorentz diffusion operator (i.e. it causes diffusion in A space) 
and 

is the gyroaveraged energy diffusion operator. 
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5.3. The Collision Operator 



UL[h\= udFo Jo {a)vii — jr— ^ — + Ji {a)v± — jr— ^ — (5.6) 

and 

/ , , f d^v Auvw Jr,( a) h , , { S'v IS.vv\J\((i]h\ 
V j d'^v Az/tiy Fo j d-^v Az/f ^ Fq / 

are the gyroaveraged momentum conserving corrections to the Lorentz and en- 
ergy diffusion operators, and 

E\h\ = UEV^Joia)Fo -^ (5.8) 

J d-^v vev Fq 

is the gyroaveraged energy conserving correction. In addition, the electron colli- 
sion operator has a piece corresponding to electron-ion collisions: 

CiM - ^S(i| (1 - f - 5 (1 + e) k + Mi„(a.)F„.)^ (5.9) 



In these equations ^ = v\\/v\s the pitch angle, Jq and Ji are Bessel functions of the 
first kind, and u\\ [hi] is the perturbed parallel ion flow velocity. Definitions of the 



velocity-dependent collision frequencies vor Az/, v\\, and ve are given in Ref. IZOy. 

The local conservation properties of the collision operator before being gy- 
roaveraged do not survive the gyroaverage, which introduces nonlocal finite Lar- 
mor radius effects. However, if the Bessel functions are expanded in perpendicu- 
lar wavenumber, the conservation properties are assured provided the first three 
moments of Cqj^\Ji\, that is, Ccxlh] with k^p set to in the Bessel functions, van- 
ish: 

/1\ 



d'^v 



,2 



C%K[h] = 0. (5.10) 



V^7 

5.3.2 Numerical Implementation 

The collision operator is implemented implicitly via two applications of Go- 
dunov splitting: starting with the result h"-^^ of advancing the collisionless gy- 
rokinetic equation (that is, the k"-^^ of Sections 14. 1 1 and |42)) , which we redefine as 
h*, first the pitch angle scattering terms are applied: 
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5.3. The Collision Operator 



h** - h* 
At 



L[h**] + UL[h* 



and then the energy diffusion terms are applied: 

- h* 



At 



D[h 



n+l] 



Un[h''^'] + E[h 



n+ll 



(5.11) 



(5.12) 



This splitting technique of treating the collision operator separately to the coUi- 
sionless gyrokinetic equation is accurate to first order in At. The combined oper- 
ation can be expressed as the inversion and multiplication of two matrices: 



h 



If* = [1 - At {L + UlT^ h* , 
= [l-At{D + UD + E)]'^ h**. 



(5.13) 
(5.14) 



The matrices L and D are kept tridiagonal through the use of a compact stencil 
for the velocity derivatives, but the matrices corresponding to the conservation 
operators are dense. Application of the combined matrices [1 — At{L + Ul)]^'^ 
and [1 — At {D + Ud + E)]^'^ can, however, be speeded up using the Sherman- 
Morrison formula, which gives x in the matrix equation Mx = b, provided M 
can be written in the form M = A + v, where where ® is the tensor product. 
We identify A with L and D, and all the conservation operators can be written 
in the form u v. And so finally, h** and h^"^^ can be calculated through two 
applications of the following formula: 



x = y2 



V2 ■ y2 

1 + V2- Z2 



Z2, 



(5.15) 



where 



2/2 = yo 



■22 — ^0 — 



-"0 ■ 2/0 


So - 


vi-yo 


1 


+ Vo ■ So. 


l + Vi- Wq_ 




Vq ■ Zo 


So - 


Vl ■ Zq 


1 


+ ^'o ■ So. 





■Wo, 



Wq. 



(5.16) 
(5.17) 



and all quantities in these formulae are defined for each application in Table 

15.11 Before the application of the Lorentz operator, the distribution function (at a 

given point in {k^, ky, 9) space) is rearranged in memory so that the A values for a 

given energy are next to each other in memory. Conversely, before the application 

of the energy operator, the distribution function is rearranged in memory so that 

the energy values for a given A are next to each other in memory. This increases 

the speed of application of the formulas. With the exception of ym, all quantities 
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Variable 


Lorentz Op. Energy Op. 


A 


1 - At (L + f/i) 1 - At {D + Ud + E) 


X 


h** 


b 


h* h** 


Ao 


1 - AtL 1 - AtD 




h'DV±Ji —Ah'V±Ji 




udV\\Jo -Az/t;||Jo 


V2 


f|l (electrons) vev'^Jq 




(ions) 


Uo 


-Atvo/du -Atvo/du 


Ui 


—Atvi 1 du —Atvi 1 du 


U2 


—Atujjvydq (electrons) —Atv2/dq 




U (10ns) 


du 


J d^v z/z)f II Fo J d^f Ai/f y Fo 


dq 


9 / r\ r 7*^ 4 7-1 

<e/2 fdhuEV^Fo 




Common 


Ai 


Ao + Uq® Vq 


A2 


Ai + Ui® vi 


Vrn 






A'r^U2 


Wm 


A^ui 







Table 5.1: Sherman-Morrison variable definitions for Lorentz and energy diffu- 
sion operator equations. 



in these formulae are time independent so they only need to be calculated once 
at the beginning of each simulation. 

Finally, the discretisation of the collision operator must be done in such a way 
as to preserve its conservation properties. This can be achieved as long as certain 
conditions are satisfied, laid out in detail in Ref. IZID, most notably that the dis- 
crete differentiation operator must satisfy a discrete version of the Fundamental 
Theorem of Calculus, and a discrete version of integration by parts. This, to- 
gether with the need for a compact stencil and the fact that the velocity grids are 
unevenly spaced limits the choice of difference scheme, allowing only first order 
accuracy. The scheme used is: 



ox ox Wj \ Xj 



where wj is the integration weight associated with xj. 
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Chapter 6 

The Implementation of Flow Shear in GS2 



6.1 Introduction 



The implementation of equilibrium flow shear is achieved by an elegant trick 
which was announced at a conference |74l] . and which is described in an unpub- 
lished note |87|1 . Since the effect of this equilibrium flow shear is central to this 
work, the implementation merits a description in some detail. 



6.2 Theory 



The system of equations (|2.10|) and (|2.12|) can be written as 

d d d d 



d d 



L 



dx' dy' dt^ dt 



9 



(6.1) 



where the operator L includes, among other things, the solution of the quasineu- 
trality condition to obtain the potential from the distribution function. 

We wish to solve this equation spectrally in both the x and y directions, but 
the explicit dependence on x introduced by the flow shear prevents this. In order 
to eliminate the flow shear we write the eikonal in terms of a time varying wave 
number fc^. — 'jEkyt: 

g — gQi[(kx-'rEkyt)x+kyy] ^ 2) 

which allows us to write 



d_ 
Ft 



(6.3) 



We have eliminated the explicit dependence on x at the cost of introducing an 
explicit dependence on t in the operator L, and hence in the response matrices 



6.3. Implementation 



(|4.33|) . This would make it necessary to recalculate the response matrices at every 
time step, which is unacceptable. We therefore define k* and g as follows |87|1 : 



(6.4) 



and 



g {k^, ky, t) 



g {kx, ky, t) 

g {k^ + ^Ekyt, ky,t) 



(6.5) 
(6.6) 



The equation is now 



dg_ 
dt 



k* k - 



9- 



(6.7) 



The explicit time dependence has been eliminated from the operator L, and the 
whole equation is now written in terms of k*, except for the fact that the deriva- 
tive on the left-hand side is at constant k^. The evaluation of this derivative is 
discussed in the next section. 

6.3 Implementation 

6.3.1 Discretisation 

We now discretise in time so that 



dg_ 
dt 



~n+l _ ~n 

At 



,kT + {l-re)kT^\ky,^^ 



[rerkx*,^ + {l-re)~glt\^.]-{(y-^) 



where without loss of generality we now let = 0. At this point we note that 
because in general A/c* is not equal to 'Je y< ky x At, L""*"^ will not be evaluated at 
^*,n+i |_,^|. ^Yie nearest discrete k* available. The difference between this discrete 
value and the actual value of k* is stored in the variable kxshif t. Rearranging 
(|6.8[l (and ignoring all other L dependences), we may write 



^,V„ = (1-AtL [kT^']) ~gi:.Uu 



(6.9) 
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and solving for (7"^^, we obtain the equation 



(6.10) 



6.3.2 Implementation 

Comparing (|6.10[) to the equivalent equation with no velocity shear: 



9-,:'={l-At"^'/'L[K]y'g]:^, 



(6.11) 



we note that the only difference is that in equation (|6.10|) , ^" is evaluated at a 
different radial wavenumber to L and ^"^^ Thus, in order to get GS2 to solve 
(|6.10|) rather than (|6.11)) , all that is necessary is an assignment step prior to the 



evaluation of (|6.11|) : 



(6.12) 



At the risk of anthropomorphising a computer code, GS2 is 'tricked' into solving 
(|6.10|) instead of (|6.11|) , because where it is expecting ^fcj.*(fn+i) on the right hand 
side, ie. ^" evaluated at the same radial wavenumber as L and g^^^ as in (|6.11|) , it 
actually gets glx*{t^) (|6.10|l . 



6.3.3 Boundary Conditions 

The effect of the assignment step (|6.12|) is to shift all the grids of the distribu- 
tion function and the fields by one grid point in the direction. At one end of 
the kx grid, the values of the grid are lost. At the other end, we choose a zero in- 
coming boundary condition. This has important implications for both linear and 
nonlinear simulations. Linear simulations should only be run with one non-zero 
value of /c*, and can only be run until this value falls off the gricil]- Non-linear 
simulations should include high enough radial modes that the perturbations are 
damped by finite Larmour radius effects before they are swept off the grid. 



^ Linear simulations with finite magnetic shear should, in fact, use the separate single mode 
implementation of flow shear which is not discussed here. 
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Part III 

A Bifurcation to a Reduced Transport 

State 



Chapter 7 
Introduction 



Chapter sUmi] are largely taken from Ref. itggi? . 



As described in Section |1.3[ a recent paper [45] which investigated the effect 
of flow shears several times larger than the linear growth rate of the ITG instabil- 
ity (one of the principal drivers of turbulence in fusion devices), found that for 
relatively low temperature gradients there is a range of flow shear values where 
the turbulence is completely quenched. The question remains how large shears in 
the equilibrium toroidal flow can be achieved, given finite sources of angular mo- 



mentum to drive the flow. In Ref. [45(1 . the turbulent flux of angular momentum 
was calculated and found to be large at moderate flow shears; if it is too large, in- 
creasing the input of angular momentum will not generate strong enough shear 
in the flow. 



Ref. |45I1 considered the standard Cyclone Base Case [|25] (see Section |3.8.4|) , 
with the normalised inverse magnetic gradient scale length s = 0.8 and a range 
of values of the temperature gradient and the velocity shear. In Part Hill which is 
published as Ref. |88,] , we consider a similar regime to that described in Ref. [|45l|] . 
with the difference that here the magnetic shear s = 0. The results reported 
in Part HII] show that the range of flow velocity gradients and ion temperature 
gradients where the turbulence is suppressed is much larger than in the case of 
finite magnetic shear. We demonstrate the existence of a transport bifurcation, 
in which a positive feedback between the increase in the flow gradient and the 
suppression of the turbulence provides the mechanism for a jump from low to 
high flow gradients, with a simultaneous jump in the temperature gradientjl] We 
show that this transition results in a state where the transport of heat is nearly 

^ Such a transition is in principle also possible with s — 0.8, but in a much smaller region of 
parameter space L89.1 



7.1. Model 



neoclassical, whilst the transport of momentum remains largely turbulent. 

In a fusion device, the quantities that can be controlled are the input of heat 
and toroidal angular momentum, which in steady state are proportional to their 
outgoing fluxes. By varying these quantities, it is possible to vary the equilibrium 
gradients. In contrast, in the local gyrokinetic simulations used in this study 
the input parameters are the gradients, and the fluxes are calculated from the 
output. In order to demonstrate the existence of a bifurcation in the gradients at 
fixed values of the input fluxes, we use local gyrokinetic simulations to map out 
the dependence of the turbulent fluxes on the temperature and flow gradients, 
and then we add the neoclassical fluxes and invert this numerically determined 
dependence to find the gradients as functions of the total fluxes. We first do this 
by straightforward interpolation in the parameter space, then propose a simple 
parameterisation of the fluxes that fits the data and can be used to predict in 
which parameter regimes bifurcations can occur. 

The rest of Part|lll]is organised as follows. In Chapter [81, we report a numer- 
ical scan in two parameters: the flow shear and the ion temperature gradient, 
calculating the turbulent heat and momentum fluxes over wide intervals of these 
parameters. In Chapter|9l, we interpolate our results and determine how the tem- 
perature and flow gradients depend on the heat and momentum fluxes: a trans- 
port bifurcation is obtained as a result. In Chapter [lOl we construct a parametric 
model of the transport. This allows us, in Chapter [TTl to study the effect of the 
transport bifurcation on the temperature gradient and investigate the range of 
heat and momentum flux values for which we expect transitions to occur. 



7.1 Model 



In Part [Till we use the gyrokinetic equation (|3.45|) to model the turbulence 



The equation is solved by the simulation code GS2, as described in PartHIl We 
vary the perpendicular flow shear 7^;, and k = Rq/Lt- Other parameters are kept 
fixed: 

^ = 2.2, go = 1.4, ^ = 0.18. (7.1) 

The density gradient, the magnetic safety factor and the aspect ratio are chosen to 
conform to the Cyclone Base Case (see Section |3.8.4|) , while, as we explained 
above, the magnetic shear is set to 0. Note that the value of go effectively controls 
the strength of the PVG drive (the term proportional to qo'jE on the right hand 
side of (|3.45|) ) for a given velocity shear 7^. The relatively low value of go in the 



Cyclone Base Case reduces the destabilising effect of the PVG — compared, for 
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example, to the Waltz Standard Case where go = 2.0 |42|1 . Collisions are included 
by means of the model collision operator described in Section 15.31 The numer- 
ical ion-ion collision frequency vj^ = O.Olvthi/Ro- All fluxes will be reported in 
dimensionless units by normalising them to the gyro-Bohm values 

The resolution of the majority of simulations was 64 x 32 x 14 x 24 x 8 x 2 
— the number of gridpoints in the k^Q and ky spectral directions, in the z spa- 
tial direction, the average number of pitch angles (i.e. the number of values of 
Hi for a given e,, which varies with the parallel coordinate), the number of grid- 
points in energy space Si and the two values of a. This grid provided sufficient 
scale separation in both spatial directions perpendicular to the magnetic field for 
calculating the turbulent transport, and sufficient resolution in velocity space to 
calculate the velocity integrals in Maxwell's equations with the required accuracy. 
A discussion of the parallel resolution is given in Chapter [ 
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Chapter 8 
Turhulent Fluxes 



8.1 Heat Flux 

Using simulations in the manner described in PartHIl the heat and momentum 
fluxes were calculated for k values in the range 4 < k < 13, vs flow shear values 
in the range < 7^; < 2. 

Fig. IS.lf b) shows the heat flux vs the ITG in the absence of flow shear. When 

= 0, the critical temperature gradient above which there is ITG-driven turbu- 
lence is K = 4.4. When k is increased above this threshold, the heat flux increases 
rapidly up to more than a hundred times the gyro-Bohm value. 

Fig. 18.21 shows the turbulent heat and momentum fluxes vs the flow shear for 
different ITG values. As the flow shear is increased from 0, the heat flux initially 
responds weakly, either increasing or decreasing slightly. As the flow shear 7^; 
approaches 1, the heat flux dips sharply for all values of the ITG. For moderate 
temperature gradients {k < 11), the turbulence is suppressed altogether and the 
heat flux drops to zero. The suppression of turbulence also happens at finite 
magnetic shears |45], but for a significantly narrower range of 



For larger 7^;, the heat flux starts to rise again. This increase is due to the 
PVG — we have verified that this effect disappears if the term proportional to 
qo^E on the right hand side of eq (|3.45[) is artificially set to 0. This revival of 



the turbulent transport at large shears was also observed in |45|1 . but was absent 
from the quasilinear study conducted in |4^, which showed the turbulence being 
completely suppressed above a sufficiently large toroidal shear. 



^ This is partly owing to the fact that at finite magnetic shear growing linear eigenmodes exist 
for non-zero flow shear, whereas at zero magnetic shear there are no such eigenmodes except 
when the flow shear is also zero; this is discussed further in PartllVI 



8.2. Momentum Flux 
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K 



Figure 8.1: ITG instability and turbulence at zero flow shear: (a) linear growth 
rate; (b) turbulent heat flux vs the normalised temperature gradient k = Rq/Lt- 
The case of k = 11, further explored in Fig. 112.11 is marked by *. 



8.2 Momentum Flux 

The momentum flux is zero when 7^; = (as might be expected in an up- 



down symmetric case, |l90t|91l]). As 7^; increases, so does the momentum flux — 
a result of the turbulent viscosity. However, when the flow shear reaches values 
at which it starts to suppress turbulence significantly, the trend is reversed and a 
remarkable situation arises in which increasing flow shear reduces the transport 
of momentum. This behaviour persists until reaches even larger values, and 
the turbulence is reignited by the PVG drive. The positive correlation between the 
velocity shear and the momentum flux is then reestablished. It is the existence of 
the window of suppressed momentum transport at moderate je and k that will 
enable the transport bifurcation analysed in Chapter |9l 
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Figure 8.2: Turbulent heat (a) and toroidal angular momentum (b) fluxes (nor- 
malised to gyro-Bohm values) as functions of flow shear for different values of 
the ion temperature gradient k; (c) turbulent Prandtl number vs flow shear (for 
cases where the heat flux is non-zero). 
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8.3 Turbulent Prandtl Number 

There is a clear correlation between the heat and the momentum flux, which 
is best quantified in terms of the turbulent Prandtl number 



^ = n^/ngg K rp ^^ ^^ 

where the turbulent viscosity and the turbulent heat diffusivity are 



T't = 7—^, (8-2) 



QgB K,2y/2Ro 

respectively. 

The Prandtl number obtained from our simulations is plotted in Fig. I8.2r c). 
There is a similar basic trend for all values of k: Fit rises from approximately 
0.5 when 7^; ~ peaks at 7s ~ 1 and then drops to approximately 0.6 when 
7e ~ 2. For those intermediate values of 7^; where the turbulent fluxes are re- 
duced or tend to 0, Pit rises sharply, reaching just above 0.7 for low values of k. 
The location of this sharp rise varies with temperature gradient similarly to the 
location where the turbulence is suppressed. 

Although the Prandtl number does vary with both 'Je and k, this dependence 
is relatively weak compared to the individual dependence of ut and Xt on these 
parameters. Thus, approximating Pr^ ~ 0.55 everywhere keeps it within approxi- 
mately 25% of the true value for the entire range of 7^;. This will prove convenient 
in constructing a model of turbulent transport presented in Chapter [TOl 



^ The low value of Pr^ for small je has been observed before and is sometimes referred to as 
a shear pinch [92, 93]. It occurs because, at low je, the perpendicular flow shear can give rise to a 
contribution to the viscous stress that has the opposite sign to that of ITG-driven turbulence with 
zero flow shear, reducing the overall diffusive transport. 
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Chapter 9 
Transport Bifurcation 



9.1 Possibility of Bifurcation 

In Chapter |8l, we demonstrated the existence of a wide interval in 7^; in which 
a sheared toroidal equilibrium flow completely suppresses turbulent transport. If 
such a suppression could be achieved experimentally in a tokamak, confinement 
of energy would be dramatically improved. Unfortunately, while it is possible 
to specify the flow shear in numerical simulations, in experiment it can only be 
varied indirectly by applying torque, and the effect of that torque is strongly de- 
pendent on how quickly the angular momentum escapes from the plasma. If 
only a limited amount of torque can be injected and the momentum flux is too 
large, it could be impossible to achieve flow shears that are large enough to sup- 
press the turbulence. Fig |8.2| fb), however, suggests an intriguing possibility. For 
all temperature gradients, the momentum flux at first increases, reaches a maxi- 
mum and then decreases before increasing again at higher flow shears. There is, 
therefore, a parameter window in which increasing the flow shear decreases the 
transport of momentum. If this were to happen, the momentum would build up, 
increasing the flow shear, which would further decrease the transport of momen- 
tum. Such a positive feedback could lead to a bifurcation and so it might seem 
that high values of flow shear could be reached without excessive input of mo- 
mentum, a possibility which was discussed in [|94l]F However, Fig |8.2| fb) shows 
the momentum flux at constant ion temperature gradient: the temperature gra- 



^ Note that although a similar mechanism was discussed in ll94ll . that study considered flow 
shears up to a maximum of approximately 0.08cs/a, where Cs = \j2Tf,lmi is the sound speed and 
a is the minor radius of the plasma (in this study ~ vthi and a = Rq). Thus they were consider- 
ing a maximum in the momentum flux which occurs at low flow shear and finite magnetic shear, 
and which is discussed in more detail in jisll . As can be seen from Ref. (i^ , the corresponding 
jump in flow shear is much smaller, and the range of k values at which such a maximum in the 
momentum flux exists is much narrower, than in the present case (Fig. l8.2r b)). 



9.2. Inverting the Problem 



dient would not, in fact, stay constant during this process. Indeed, as soon as 
flow shear began to suppress the turbulent transport of momentum, it would 
also suppress the turbulent transport of heat, causing an increase in the temper- 
ature gradient as well as in the flow gradient. This increase in the temperature 
gradient would restore the turbulence to its former levels. Nevertheless, we will 
show in this section that when neoclassical transport is also taken into account, a 
bifurcation is possible. 



9.2 Inverting the Problem 

What can actually be controlled in a steady-state experimental situation is the 
flow of heat and momentum through a particular surface. This means that we 
must switch from using 'Je and k as independent parameters to using the total 
fluxes of heat and momentum, Q and 11 respectively. 

Let us consider an experimental set-up where the heat and the momentum 
are injected by beams of neutral particles. We assume that the energy and mo- 
mentum from those neutral beams are deposited uniformly across the torus. We 
further assume that the beams are tangential to the magnetic axis and we ignore 
corrections of order the aspect ratio of the device. Then: 



Q _ 1 VfN,myj2 _ 2V2i?oPNBi 



QgB ~ QgB 47r2roi?o An'^roTiiTiVthip' 

JI 1 VfNbmjVb 

UgB UgB 47rVo 

Q/QgB \/2Vb v2-EnBI 



(9.2) 
(9.3) 



where Nb is the number of neutral beam particles injected per unit time, Vb is the 
beam particle velocity, Pnbi = Nbniivl /2 is the beam power, -Enbi = raivl /2 is the 
beam particle energy and Vj is the volume fraction of the plasma enclosed by the 
flux surface. Thus, the total heat flux Q is determined by the beam power Pnbi 
(|9.1|) , whereas the ratio of the total momentum angular flux to the total heat flux, 
Ti/Q, is determined by the beam particle energy Pnbi (|9.3[) . We want to know 
whether, by varying Q and n/Q, it is possible to reach, or to trigger a transition 
to, a high-flow-shear regime where the turbulent transport is suppressed. Thus, 
it is necessary to convert the dependence of Q and 11 on k and 7^, determined 
from local simulations, to a dependence of k and 7^ on Q and li/Q. 
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9.3 Interpolation 

The gyrokinetic code GS2 gives the fluxes as a function of k and 7^;. Given the 
expense of the nonlinear simulations necessary to do this, it is computationally 
challenging to use GS2 as a root finder to invert the problem and find the gradi- 



ents as a function of the fluxes |85l,l95|1. Instead we interpolate within the set of 
data points described in Chapter |8] (as well as additional simulations in parame- 
ter regions of particular interest) to obtain the fluxes as continuous functions of 
the gradients; these functions can then be inverted numerically to give k and 7^; 
as functions of the fluxes. 

Interpolation in a multidimensional parameter space is a nontrivial operation. 



A standard technique is to use radial basis functions \96\. which weigh each data 
point by its distance in parameter space from the point of interest (after the pa- 
rameter space has been normalised to ensure that variation occurs on the same 
scale in each coordinate). There are many choices of the function, or kernel, which 
is used to calculate the relative importance of each data point. Here we choose 
a linear kernel (equivalent to linear interpolation in the case of only two points) 



196(1 . Using this, the values of the fluxes at a point p = (7E, can be calculated 



as follows: 



Q,(p) = ^^f (9.4) 
j 

nt(p) = 5^<|p-Pj|, (9.5) 

j 

where pj are the input gradients for the nonlinear simulation labelled j, and the 
weights and are calculated so as to satisfy for all j: 

gi(Pj) = Qij,n,(pj) = n,j, (9.6) 

where Qtj and liti are the values of the fluxes obtained from simulation]. 

9.4 Neoclassical Transport 

If the turbulent transport is successfully suppressed, neoclassical (coUisional) 
transport becomes important. The total fluxes are the sum of the neoclassical and 
turbulent contributions: 



Q = + Qn, n = + n 



(9.7) 
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9.5. Bifurcation 



where the neoclassical fluxes are 



Qn 2\/2Rq 

= Xnfi 2~' (9.8) 



= VnlE 2' (9-9) 



l-gB 

and the neoclassical thermal diffusivity and viscosity are 



R 



3/2 



Xn ^ 0.66 ( ^1 qlp^uu, (9.10) 



z/„ ~ O.lqy.uu. (9.11) 

The formulae for the neoclassical diffusivity and viscosity in the large-aspect- 
ratio, concentric-circular-flux-surface, banana {vu <^ vthi/Rq) regime, were taken 



from [ISy. In |15|1 , the ion-ion collision frequency is defined as 

~ ^372 172 



(9.12) 



(where In A is the Coulomb logarithm). Thus, the neoclassical transport is a func- 
tion of the temperature and density, and scales with these quantities in a different 
way to the turbulent transport. To determine the ratio between the neoclassical 
and turbulent transport it is necessary to choose specific values for the tempera- 



ture and density. Here, as in [46(1 . we take uu = z/^^ /2. Later, when we consider 



the case of specific tokamaks, we will use typical values of Ui and Tj to estimate 



Using the results (|9.10|) and (|9.11|) the neoclassical Prandtl number can be cal- 



culated from the parameters given in (|7.1[) and is found to be 



Pr„ = — ~ 0.01 < Pif (9.13) 

Thus the neoclassical Prandtl number is much smaller than the turbulent Prandtl 
number (see Fig. |8.2| fc)). In other words, the neoclassical transport of momentum 
is much smaller than the neoclassical transport of heat, in contrast to turbulent 
transport, for which the fluxes of momentum and heat are comparable. While 
the formulae (|9.11|) and (|9.10|) are approximate, we emphasise that the qualitative 



result of the following section is not affected by small changes in the values of 
Xn and z/„, provided that the property Pr„ <^ Fit continues to hold. In particular, 
this means that this qualitative result is not affected by small changes in the value 

of Uii. 
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Figure 9.1: Momentum flux divided by the heat flux (each normalised by the 
corresponding gyro-Bohm estimate) vs the flow gradient at a constant value of 
the heat flux Q = 2.QQgB plotted using (a) interpolation from the data explained 
in Section |93l and (b) the parameterisation from Chapter [TOl Also shown are the 
neoclassical contribution to the momentum flux (dotted line) and the momentum 
flux at constant heat flux without the neoclassical contribution (dashed line). 



9.5 Bifurcation 

In Fig. I9.ir a), the momentum flux is once more plotted against 7^;, but this 
time keeping the heat flux constant at Q/QgB = 2.6. The local maximum in the 
momentum flux still exists. Starting at this local maximum (point A), if the torque 
on the plasma was to be increased at constant Q, the flux of momentum could 
only increase if the flow shear were to jump to a much higher value (point B) 
where the PVG instability would drive turbulent momentum flux. A bifurcation 
is manifest. 

The inclusion of the neoclassical fluxes is critical in obtaining this result. To 
illustrate this. Fig. 19.11 also shows the momentum flux at constant heat flux with- 
out the neoclassical contribution to the fluxes, i.e., with Xn = J^n = 0. If, in such 
a synthetic "purely turbulent" transport case, the flow shear is increased at con- 
stant Q, the temperature gradient increases to maintain the turbulent heat flux. 
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9.5. Bifurcation 



and, as shown in Fig. 19.11 Ti/Q increases monotonically with 7^;. The key prop- 
erty of neoclassical transport that helps change this behaviour and produce a 
bifurcation is that the neoclassical Prandtl number is much smaller than the tur- 
bulent Prandtl number. This means that as turbulent transport is suppressed, 
the system goes through a regime where the neoclassical contribution to the heat 
flux is significant, while the neoclassical contribution to the momentum flux is 
not. So as n/Q is increased at constant Q, the turbulence is suppressed, and the 
temperature gradient starts to rise, but as this happens more of the heat flux is 
transported neoclassically, and so the feedback loop breaks down: it is no longer 
necessary to increase the levels of turbulence to maintain the same Q. The same 
is not true of the momentum: the neoclassical viscosity cannot make up for the 
lack of turbulence, and so the transport of momentum peaks and then falls. 

At high flow shears the turbulent momentum flux increases rapidly once 
again; this turbulent flux is driven by the PVG, as discussed in Chapter |8l As 
a result, we observe that the bifurcation results in a turbulent state at B, with a 
significant turbulent momentum flux. This is a substantial difference from the 
situation envisioned in |40|l , where a transport bifurcation was predicted using a 
reduced quasilinear model. Their model did not (and, being a quasilinear model, 
could not) predict the existence of the PVG-driven subcritical turbulence at high 
flow shears; instead it predicted a bifurcation resulting in a non-turbulent state 
where all transport was neoclassical. Thus a full nonlinear analysis is necessary 
to describe accurately the reduced transport state produced by the bifurcation 
that we find in a turbulent plasma. 

Restoring dimensions, we find that the bifurcation results in a jump in the 
flow shear from 0.25 Vthi/Ro to ^-^7 Vthi/Ro- It is instructive to consider whether 
such values of shear appear in current devices. Taking the measured values of the 
perpendicular flow shear, Tj and Rq, in high-performance discharges in the JET 
and MAST tokamaks |52l,l97|1 , we estimate that flow shears of up to approximately 
Vthi/Ro have been recorded in both devices, and thus that values of shear of the 
order examined in this chapter have been observed. 
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Chapter 10 
Parameterised Model 



Through simple interpolation, without making any assumptions about the 
way the fluxes depend on the gradients, we have shown the existence of a bifur- 
cation in the gradients at constant fluxes. However, even interpolating one line 
of constant Q requires a very large number of data points in the region of the 
line. To produce the required data set for Fig. I9.ir a), we had to perform about 350 
nonlinear gyrokinetic simulations, each run to saturation. In order to extend our 
understanding of the bifurcation, and of the range of flux values for which similar 
bifurcations can occur, we will first consider the behaviour of the turbulent fluxes 
in further detail and then use this analysis to construct a simple parameterised 
model of the turbulent fluxes Qt and lit as functions of 7^; and k. 



10.1 Modelling Qt 

Dependence of Qt on k 

We wish to construct a parameterised model of Qt as a function of 7^; and 
K. In Chapter |8] we described the dependence of Qt on 7^;; we now consider 



the dependence of Qt on k. Both experimentally [98'] and numerically [25(1 . it is 



usually found that Qt increases very sharply with k — a property known as stiff 



transport. A recent experimental study 19811 has indicated that increasing the flow 



shear at low magnetic shear might have the effect of reducing the stiffness. Fig. 
110. II shows that the effects of flow shear on dQt/dn are in fact quite complex. Let 
us consider the cases of low flow shear, 7^; < 1, and high flow shear, 'Je > 1, 
separately. 

For the case of 'Je < 1/ shown in Fig. llO.ll fa), the threshold in k above which 
turbulence can be sustained nonlinearly increases rapidly with 'Je, as the per- 
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Figure 10.1: Turbulent heat flux vs the temperature gradient for different values 
of the flow shear, showing (a) low-shear values 7^; < 1, (b) a close up of the low 
Qt region in (a), and (c) high-shear values 7^ > 1 
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pendicular shear suppresses the ITG instability. Above this threshold there are 
broadly three regions: low, intermediate and high Qt. 

At high Qt, far above the threshold, for any value of flow shear the heat flux 
eventually asymptotes to the same dependence as it has at 7^ = 0. In other 
words, as the ITG drive becomes very strong, the effect of flow shear becomes 
negligible. 

At intermediate Qt, the heat flux rises rapidly from low values to join the 
universal, high-Qt asymptotic. As 7^; increases, because the threshold rises, the 
heat flux rises more rapidly with k above the threshold; thus at intermediate 
values of Qt, the stiffness dQt/dn increases with 7^;; this is shown in Fig. I10.2[ 

At low Qt, just above the threshold, the heat flux rises very slowly, that is, the 
stiffness dQt/dn is very low (see Fig. I10.2[) . This low Qt region is shown in detail 
in Fig. 110. ir b). Such is the sharpness of the transition from the low-stiffness low- 
Qt region to the high-stiffness tntermediate-Qt region that there appear to be two 
distinct thresholds. The first threshold is the transition from no turbulence to 
non-zero turbulent transport. Above the first threshold turbulence is present but 
Qt rises slowly with k (the \ow-Qt region). Above the second threshold Qt rises 
rapidly (the tntermediate-Qt region). These thresholds are plotted in Fig. |10.3[ 
The low-stiffness region only exists for < 7^; < 0.8: between 0.4 < 7^; < 0.8 the 
first threshold joins the second and the low-stiffness region disappears. 

At high flow shear, 7^; > 1, there are two principal differences to the case 
with low flow shear. Firstly, the low-stiffness region is not present, and there is 
only one threshold. Secondly, the critical temperature gradient for turbulent heat 
transport starts to decrease as the PVG drive re-enforces the ITG drive. When 

= 1-8, the PVG drive is strong enough to drive turbulence at very low k: 
the threshold drops to zero. As the threshold drops, the transport stiffness at 
intermediate values of Qt decreases, in a mirror image of the case at low flow 
shear. At high Qt, the heat flux still asymptotes to the universal 7^ = curve. 

Thus, increasing 7^; can both increase and decrease the nonlinear thresholds, 
and both increase and decrease the stiffness. The rise and fall both in the thresh- 
olds and the stiffness can also be seen in the finite-magnetic-shear simulations of 
1 4^. The principal differences between the zero-magnetic-shear case considered 
here and that case are that we have a higher value of the critical gradients for all 
7^;, and that we find a low-stiffness region at low values of — a feature that 
seems to be absent at finite magnetic shear. 



80 



10.1. Modelling Qt 



60 

50 

!£ 40 
> 30 
^ 20 
10 




1 1 1 1 1 1 


1 1 










Low Qt 

: 1 t t ^ 


1 1 





0.2 0.4 0.6 0.8 1 

IE 



1.2 1.4 1.6 1.8 



Figure 10.2: The dependence of the heat transport stiffness dQ/dn on the flow 
shear in the \ow-Qt and intermediate-Qt regions, measured using simulations 
close to both thresholds (points, see Fig. |10.1|) and using the parameterised model 
of Section lTO.il (lines). 
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Figure 10.3: The first and second critical thresholds, k^x and Kc2, measured using 
simulations close to both thresholds (points, see Fig. 110. 1|) , and using the param- 
eterised model of Section [1 . 1 1 (lines) . 



Parameterisation of Qt 

In Ref . 1890 ^ simple model was used to characterise the behaviour of the tur- 
bulent heat flux, and describes the qualitative behaviour of the heat flux reason- 
ably well. Here, we describe a more complex model that reproduces most of the 
features described above and is quantitatively close to the interpolated fluxes. We 
consider only the low and intermediate Qt regions, as the bifurcation occurs near 
the boundary between these regions. In order to describe Qt in these regions, we 
have to parameterise the two thresholds and the transport stiffness dQt/dn. 

We parameterise the first and second critical thresholds as linear and 
quadratic functions of 75; respectively: 
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i^ci = 1^0 + ailE, (10.1) 
Kc2 = /to + a27E + "371, (10.2) 

where kq = 4.4 is the nonlinear threshold for turbulence at 7^; = and the param- 
eters tti = 10.1, a2 = 17 A and = —17.0 are chosen to fit the data. The modelled 
thresholds are plotted next to the measured thresholds in Fig. 110.31 As with the 
observed thresholds, the first threshold joins the second between 0.4 < 7^; < 0.8. 

Next we parameterise the transport stiffness dQt/dn. The measured values 
of dQt/dn in the first and second regions are shown in Fig. 110.21 Between the 
first and second thresholds (the low Qt region), d{Qt/QgB)/dK is modelled as a 
constant, 04. In the intermediate Qt region, remembering the observation that 
dQt/dn broadly rises and falls with the nonlinear thresholds, we allow the stiff- 
ness to depend on Kc2- Thus, 

1 dQt J "4 K^i<K< K,2, 



QgB dn I as + aQKc2 ^ > /«c2- 

The parameters are = 1.5, = 3.0 and = 8.0. The model of dQt/dn is 
shown in Fig. 110.21 

Thus Qt is parameterised as a piecewise linear function of k. It is zero below 
the first threshold, has gradient above the first threshold and gradient as + 
aQKc2 above the second: 



Qt 



K < Kcl, 

^ a4 (k - Kcl) Kcl < K < Kc2, (10-4) 



QgB 

The model is compared with the data points in Fig. |10.4[ To reproduce the bifur- 
cation of Section 19.51 in a quantitatively correct manner, it is in fact sufficient to 
represent accurately the two thresholds and the low stiffness region. However, 
by also parameterising the variation of dQ / dn, we have provided a model which 
uses six parameters to describe the (complicated) behaviour of the heat flux over 
a wider parameter regime with reasonable accuracy. 



10.2 Modelling 

While the turbulent Prandtl number does vary with 7^ and k, as discussed in 
Section 18.31 this variation is relatively weak. Thus we model Ut by assuming a 
constant turbulent Prandtl number: 
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Figure 10.4: The modelled heat flux (lines) shown along with the simulated data 
points for (a) low flow shear and (b) high flow shear. Legends as in Figures flO-lT a) 
and Sane). 
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(10.5) 



The same choice was made in Ref . |89[] 



10.3 The Modelled Bifurcation 

In Fig. |9.1| fb) the model is used to replot the angular momentum flux at con- 
stant Q/QgB = 2.6. It has the same shape as the interpolated curve in Fig. 19.1 I f a), 
and is quantitatively very close to it at low flow shear (7^ < 1.0). The agree- 
ment at higher flow shears is less good, which reflects the fact that the second 
critical gradient is not really a quadratic. Nonetheless, we consider this degree of 
precision adequate. The model is used in Chapter [TT1 to explore further proper- 
ties of the transition, without the need for extremely large numbers of additional 
simulations. 
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Chapter 11 
The Reduced Transport State 



In Chapter m a transition was described leading to a reduced transport state. 
Here we use the parametric model that was designed in Chapter [10] to describe 
the properties of this state. 



11.1 Heat Flux at Constant n /Q 



Fig. 19.11 showed the effect of varying H/Q whilst keeping Q constant. In con- 
trast. Fig. II 1 . 1 1 demonstrates the effect of varying Q whilst keeping H/Q constant. 
At high Q and low temperature gradient (marked (I) in Fig. Ill.ll fb)), increasing 
the heat flux has the effect of increasing the temperature gradient as might be 
expected: this is the ordinary turbulent regime. However, as Q/QgB is lowered 
below ~ 5 there arises a counterintuitive situation where decreasing Q has the 
effect of raising the temperature gradient. This is of course the combined effect 
of the neoclassical Prandtl number and the flow shear as described in Chapter 



|9l, and also discussed in detail in In this region (marked (II) in Fig. Ill.lf b)), 
Qt becomes comparable to Qn, but as the heat flux is lowered, the system cannot 
drop to a neoclassical state because of the smallness of the neoclassical transport 
of momentum. Thus the constancy of the ratio li/Q (i.e., the large finite flux of 
momentum) keeps the system in a turbulent state; the temperature gradient and 
the flow shear rise, and the curve representing Q vs k becomes flatter and curls 
up in order to maintain its distance from the neoclassical line. At large temper- 
ature gradients and flow shears (region (III) in Fig. IlLlT b)), the momentum flux 
increases rapidly due to the PVG drive and the high flow gradient (see Fig l8.2l fb)) 
and so the curve rolls over and resumes its original downward trend. Eventually 
the heat flux asymptotes to the neoclassical value. 

The last part of this trend, while physically obvious, can only be seen using the 



11.1. Heat Flux at Constant U/Q 




Figure 11.1: Heat flux Q vs the temperature gradient k at a constant ratio of the 
momentum flux to the heat flux U/Q (in units of \/2Ro/vthi) plotted using (a) 
interpolation from the data (Sec 19. 3[) and (b) the parameterised model (Sec ITOl) . 
Also shown is the neoclassical contribution to the heat flux. Interpolation in (a) is 
impossible near this neoclassical line, where the contours are closely spaced and 
U/Q is multivalued, (b) also shows the maximum possible temperature gradient 
at low Q for a given U/Q (labelled "max k"). Fig. (c) plots the same curves 
as (b), showing both Q against 7^; and k, and the curves projected on the 7^;- 
K plane, illustrating the increase in the flow shear along each curve of constant 
U/Q. Points A and B in (a) correspond to points A and B in Fig. |9.ir a). 
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11.2. Temperature Gradient Jump 



modelled heat flux (Fig IlLlT b)), as it is not feasible to interpolate in the region 
near the neoclassical line where the contours of constant Ti/Q are very closely 
spaced and li/Q becomes multivalued. 



11.2 Temperature Gradient Jump 

Interpolation cannot yield the temperature gradient after the transition di- 
rectly: as explained above, the contours of constant Ti/Q become too closely 
spaced as they approach the neoclassical line in Fig. Ill.lT a). However, the tem- 
perature gradient of the final state, labelled B on both Fig. 19.11 and Fig. Ill.lf a.b), 
can be calculated indirectly by using the value of 7^ from Fig. I9.ir a) and rear- 
ranging (|8.1|) and (|9.7[) to give k as a function of 7^;, Pif, li/IigB and Q/Qqb'- 



Ti/TigB + {4:Rlqo/vthiropf) 7^; (x„ - «^n) ' 
Taking 7^ = 1.17, Pr^ = 0.55 yields a temperature gradient at point B of k = 9.7. 
The temperature gradient at point A (the other side of the jump) is 7.4. Alter- 
natively, using our model, the values of the temperature gradient can be read 
from Fig. Ill.ll fb) which results in an identical jump of k = 7.4 to k = 9.7. If we 
compare with the case of zero momentum input and zero flow shear (point C on 
Fig. Ill.ll fa)), we find that flow shear has enabled a total increase in the tempera- 
ture gradient of a factor 9.7/4.5 = 2.2. If this were reproduced across the whole 
device, it could increase the core temperature bva factor of around 10, but this 



would require a more detailed ID model (see ISSD), rather than the OD model 
presented here. Note that while the jump in temperature gradient caused by the 
bifurcation (the transition from A to B) is a striking feature, a comparable con- 
tribution to the increase of the temperature gradient arises from the incremental 
suppression of the turbulent transport by the velocity shear (the difference be- 
tween A and C). 



11.3 Neoclassical Heat Flux, Turbulent Momentum 
Flux 

It was noted earlier that the reduced transport state produced by the bifur- 
cation is still turbulent, with a momentum flux far greater than its neoclassical 
value (Fig. 19. 1|) . In Fig. Ill.ll fb), it can also be seen that the reduced transport 
state does not lie exactly on the neoclassical line; it does however lie very close 
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11.4. Bifurcation by lowering Q / QgB 



to it, which implies that the turbulent heat flux is much smaller than neoclassi- 
cal. Thus, the bifurcation takes the system into a state where the heat is mostly 
transported neoclassically, but the angular momentum is mostly transported by turbu- 
lence. The small ratio of Qt to Qn reflects the fact that the turbulence levels in 
the reduced transport state are small. The dominance of 11^ over n„ given the 
low levels of turbulence reflects the fact that the flow gradient is large and the 
neoclassical momentum transport is very low compared to the neoclassical heat 
transport (Pr„ ^1). 



11.4 Bifurcation by lowering Q / QgB 

Starting from point A, if Q were to be reduced at constant U/Q, the system 
would again be forced to jump to point B. In effect, what would happen is that the 
heat flux would become principally neoclassical; the momentum flux would drop 
because Pr„ <^ Pr^, and a bifurcation would ensue in the manner described in 
Section 1931 Thus, a decrease in the input heat could lead to a higher temperature 
gradient. We note, however, that Q is normalised to the gyro-Bohm value: QgB = 
niTiVthipff^v^Rl. Thus, decreasing Q/Qgs could correspond to decreasing the 
deposition of heat, but it could also correspond to an increase in temperature or 



density ||99t]. 



11.5 Optimum Q 

Fig. Ill.ir b) shows that for every value of U/Q, there is an optimum Q that 
gives a maximum temperature gradient (the dotted line at k ~ 11 in Fig. IlLlT b)). 
The maximum temperature gradient increases with U/Q, but only slowly. The 
optimum value is the value of Q such that: 



dn 
dQ 



= (11.2) 

n/Q 

and is plotted as a function of n/Q in Fig. 111.21 The maximum temperature gra- 
dient is effectively determined by the maximum nonlinear critical temperature 
gradient Kc (see Fig. I10.3|) . In Chapter [14] we will show how this maximum of Kc 



can be dramatically increased. 



87 



11.6. Transition Region 



Table 11.1: Typical plasma properties from the ITER Profile Database |100(1 . The 
symbol a denotes the minor radius of the device. The temperature was calculated 
as the mean of the core (TIO) and edge (TI95) temperatures. 





T, (eV) 


rii (m ^) 


R{m) 


Bt{T) 


a (m) 


DIII-D 


2.7e+03 


9.0e+19 


1.7e+00 


1.7e+00 


6.1e-01 


JET 


2.6e+03 


6.3e+19 


2.9e+00 


2.4e+00 


9.3e-01 


MAST 


0.6e+03 


3.9e+19 


8.0e-01 


5.3e-01 


5.6e-01 



11.6 Transition Region 

Finally, we will show that there is in fact only a limited range of both Q and 
H/Q where bifurcations can happen in the way described above. To illustrate, we 
observe that no transition can occur along the contour U/Q = 0.12\/2Ro /vthi in 
Fig. lll.ir b): if n/Q is kept constant at this value, decreasing Q/QgB from greater 
than its optimum value of 4.0 increases the temperature gradient smoothly up to 
its maximum value. The existence of a bounded region where such transitions 
can happen is studied in detail by the authors of Ref. ||89tl , who calculate, using 
a simple transport model, the region in which transitions occur, and derive a 



criterion necessary for their existence. We apply the analysis of Ref. |89|1 using the 
model for the turbulent fluxes given in (|10.4|) and (|10.5|) , to determine the range 
of values of Q and H/Q for which bifurcations of the form we have described can 
occur. In essence, bifurcations can only occur when there exist multiple values 
of jE and K that give rise to the same values of Q and H/Q. From Fig. Ill.lf b), it 
is clear that this is only possible for values oiU/Q where there is a minimum in 
the curve of constant H/Q, and for values of Q which lie between this minimum 
and either a maximum in the curve or the point where the curve intercepts the 
neoclassical line. Thus the region in which transitions are possible is bounded by 
the curve 



dQ 
Ok 



= 0, (11.3) 

U/Q 

and by the neoclassical transport curve. This region is plotted in Fig. 111.21 

In order to give a clearer meaning to this diagram, we use equations (|9. 11 - 19. 3[) 
and typical properties of plasma devices taken from the ITER Profile Database 



llOOy (and listed in Table lTLTj) to replot the region in which transitions can happen 
in terms of the input beam power and beam particle energy, Pnbi and -Enbi- To 
generate these plots, we also calculate the collision frequency ua self -consistently 



using (|9.12|) , to take account of the varying strength of the neoclassical transport 
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11.7. Summary 



in each device. The transition regions for each device are displayed in Fig. 111.31 
along with a point indicating the typical values of Pnbi and -Enbi in each device. 
Fig. lll.3l shows that in which transitions can happen lie within an order of magni- 
tude of the typical values. It should be stressed that with the simplified magnetic 
geometry (Chapter (T)) and model of neoclassical transport (Section |9.4[) used in 
this study, and with the many assumptions about the way the energy and mo- 
mentum are deposited in the plasma (Section I9.2[l , closer agreement could not 
be expected. In particular, the assumptions of Section 19.21 are likely to overesti- 
mate the applied torque and hence overestimate the beam energy needed for a 
transition. 

4 




0.06 0.08 0.1 

^/Q/(i?oMh^) 

Figure 11.2: The region in which bifurcations can occur, calculated using the anal- 
ysis of Ref . |89|1 and the parameterisation of the fluxes described in Chapter [lOl 

The 



The cross indicates the location of the bifurcation described in Section 19. 5[ 
dashed line represents the optimum value of Q for a given value of Ti/Q, Eq. 
(fTLlb . 



11.7 Summary 

We have demonstrated the existence of a transport bifurcation to a high tem- 
perature gradient, high flow gradient state, and estimated the values of heat and 
momentum input at which we might expect such a transition to occur. In the 
reduced transport state the heat flux is nearly neoclassical but the momentum 
flux is large owing to the existence of turbulence driven by the PVG. It is this 
turbulence was effectively limits the maximum temperature gradient that can be 
reached in the bifurcation: owing to the stiffness of turbulent transport, this max- 
imum gradient {k = 9.7) is, in fact, close to the maximum nonlinear critical tem- 
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Figure 11.3: The regions (shaded) in which transitions can happen, as a function 
of the beam power and the beam particle energy, for three different devices. The 
dashed lines indicate, for each device, the value of Pnbi for a given -Enbi which 
would lead to the optimum temperature gradient, as described in Section |11.5[ 
The points indicate typical values of Pnbi and Pnbi for ea ch d evice. The projected 
Pnbi and Pnbi for MAST Upgrade were taken from Ref. flOll] . 



perature gradient for turbulence (kc ~ 11, see Fig. I10.3|) J^I In Part|IV]we consider 
this PVG-driven turbulence in detail, and show that the maximum temperature 
gradient that can be reached without destabilising it can be greatly increased. 



The difference between the two is owing to the neoclassical transport. 
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Part IV 
Subcritical Turbulence 
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Chapter 12 
Turhulence without Instability 



At zero flow shear, turbulence driven by the ion temperature gradient results 
from the presence of unstable eigenmodes in the system. The instability that is 
present in the system means that arbitrarily small perturbations to the electric 
field and the particle distribution function can grow to the size needed for them 
to interact and form turbulence. As was discussed in Section [L3l a recent paper 
1 I45I1 which investigated the effects of flow shear on turbulence driven by the ion 
temperature gradient — at a finite value of magnetic shear {s = 0.8) — showed 
that at low non-zero values of the flow shear, unstable eigenmodes were present 
in the system considered (see also the earlier results in slab geometry of Artun. 
et. al. 1480), but that for values of 7^; higher than approximately 0.8, no such 
eigenmodes existed. 

In contrast to the results of Ref. |45l] . we find in the case of zero magnetic 
shear that the growth rate is zero for all non-zero values of 7^;. With 7^; = 0, the 
ITG instability exhibits robust linear growth, with a threshold of k = 4.4 (Fig. 
I8.ir a)). However, for any finite value of the flow shear perturbations grow only 
transiently before decaying (Fig ll2.1l a)). While formally this means that 7^ — ^ 
is a singular limit, there is in fact no physical discontinuity in behavior: as 7^; — )■ 0, 
the duration of the transient growth, t^, tends to infinity (roughly as l/7£;; see Fig. 
112. If b)). Thus, there is a continuous transition at 7^; = from a transient mode 
that grows for an infinitely large time, to a growing eigenmode. 

The differences between the two cases may be explained as follows. At finite 
magnetic shear, growing eigenmodes are possible because while the perpendic- 
ular flow shear leads to an effective dependence of the radial wavenumbers on 
time (see equation (|6.2[) ), the magnetic shear makes them dependent also on the 
position along the field line. Therefore, for modes moving with a speed propor- 
tional to 7£;/s, the magnetic shear cancels the effect of the velocity shear on (t) 



(as was also noted in Refs. |39l, |49t]). When the magnetic shear becomes very 
small, or the flow shear very large, the required mode velocity becomes much 
greater than the thermal speed of the ions and the cancellation is no longer pos- 
sible because the mode cannot travel so fast In this situation, the radial 
wavenumber increases inexorably until the perturbations are damped, either by 
collisional dissipation in the fluid-like model of Ref. |49[] or by finite-Larmor- 
radius effects in the gyroktnetic model (also used in this work) of Ref. |45l] . As a 
consequence, growing linear eigenmodes only exist for non-zero magnetic shear 
and only up to a finite value of the flow shear |45l] . 
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Figure 12.1: Transiently growing linear modes at k = 11. (a) The heat flux as 
a function of time, normalised to its value at to, where to = O-OQV^Ro/vthi, so 
chosen to skip the short initial transient associated with a particular choice of 
initial condition. All modes grow transiently for 'Je > 0. (b) Duration of transient 
growth Ty from t = to to the peak value of Qt vs flow shear. 



In a standard picture of ITG turbulence with 7b = |25(1 , the turbulence is 
driven by unstable linear modes. With the exception of a narrow interval of 
temperature gradients where self-generated zonal flows suppress the turbulence 
(the Dimits shift), the presence and the amplitude of the turbulence are largely 
determined by the presence and the growth rate of those unstable modes. In 
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the present case, by contrast, we have strong levels of turbulence sustained non- 
linearly in a parameter regime where there are no linearly unstable modes, a 



phenomenon first discovered in tokamak turbul ence in R efs. |45l] and [46(1 , but 



well known in various hydrodynamical contexts |102l,ll03ll . This phenomenon is 
known as subcritical turbulence. 

Subcritical turbulence differs from standard instability-driven turbulence in 
several important ways. Firstly, because there is no linear instability, turbulence 



will not grow from initial perturbations of arbitrarily small amplitude |45l1 . Fluc- 
tuations must be initialised with sufficient amplitude (generally of the order of 
their amplitude in the saturated state) in order for turbulence to be sustained; 
thus, the absence of sustained turbulence in a numerical experiment may merely 
indicate that the initial amplitude is insufficient, not that the plasma is quiescent. 

The second difference concerns the scales at which the turbulent energy re- 
sides. In ITG-driven supercritical turbulence (7^ = 0), these scales are those 
where the linear drive injects energy — this tends to correspond to A; go-Ro ~ 1 
and kyPi relatively narrowly concentrated around values of a fraction of unity (at 



least for low values of go arid moderate temperature gradients; see Ref. |104i1 and 
Fig. I12.2l fa)). In subcritical turbulence, the preferred wavelengths appear to be 
those at which the amplification of the transient modes is strongest. In the case 
of turbulence where the PVG drive is dominant, we will show in Chapter [131 that 
the amplification is maximised along a line in Fourier space where 

Fig. 112. 2r b) shows the spectrum of strongly PVG-driven turbulence at high flow 



shear (7^; = 2.2). We see that although the result (|12.1|) is based on linear theory 
and is derived for slab geometry, it appears to describe the peak of the spectrum 
reasonably well. The spectrum extends to much higher parallel wavenumbers 
than in the standard ITG case; consequently, higher parallel resolution is neces- 
sary to resolve it\ 

Finally, faced with subcritical turbulence, we are left without an intuitive way 
of estimating its saturation level and, consequently, the turbulent fluxes. It is the 
presence of linear eigenmodes with a defined growth rate (positive or negative) 



^Studies of the effect of increasing the parallel resolution showed that at values of flow shear 
^ 1-2, the parallel resolution used for the parameter scan described in Chapter [5] (14 grid- 
points) was insufficient and led to errors in the critical temperature gradient (above which sub- 
critical turbulence could be sustained) of order 5-10%. We do not consider this error significant 
enough to merit a repeat of the parameter scan with substantially higher parallel resolution (al- 
though in principle such an exercise would be useful). 
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Figure 12.2: The spectrum of turbulent fluctuations for normal and subcritical 



turbulence: Yl 



(arbitrary units) vs ky and k\\ for (a) ■je = 0.0, k = 12 



(ITG turbulence with no flow shear) and (b) 'Je = 2.2, k = 12 (subcritical, strong 
PVG-driven turbulence) E|. Also shown in (b) is the line of maximum transient 
amplification (|12.1|) . 



which has enabled the quasilinear modelling of the heat fl ux in many situations 
without res orting to full nonlinear simulations (see Ref. |105ll for an overview 
and Refs. |106l. Il07|1 for recent work). When the growth of all modes is tran- 
sient, such models will not work in their current form. The question arises as to 
which characteristics of the linear transient growth are relevant in the resulting 
nonlinear state. 

In order to answer this question and, more importantly, the question of 
whether it is possible, without using nonlinear simulations, to predict whether 
subcritical turbulence will be present for a given set o f parameters, we present in 
Chapter [13] a detailed study — published as Ref. |108|1 — of the transient growth 
of the fluctuations. To make analytic solutions obtainable we work in a simplified 
slab geometry with straight magnetic field lines and periodic parallel boundary 
conditions which nevertheless contains the key physical effects of the ITG, the 
perpendicular flow shear and the PVG. In such a system we may Fourier trans- 
form the fluctuations in the parallel direction so that each mode has a parallel 
wavenumber k^. We show that, defining the origin of time such that A;^ = at 
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the (now time-dependent) growth rate of the fluctuations is equal at t = 
to the growth rate of an unsheared mode with the same wavenumber. As time 
increases, the growth rate of the fluctuations eventually reaches zero, at a time 
to which can we calculate in certain limits. For times greater than the mode 



To characterise the vigour of this transient growth we define the amplification 
exponent 



such that is the ratio between the size of the perturbed electric field at the 
time when it starts to decay and its size at time zero. The quantity 7(t) is the 
time-dependent growth rate. We then define A^max/ the maximal amplification 
exponent, to be maximised across all ky and k\\. We find that in general A^max 
depends on three parameters, the ratio between the parallel and perpendicular 
flow shears, q/e, the ratio between the ion temperature gradient and the perpen- 
dicular flow shear, r]s, and the ratio between the ion temperature gradient and the 
ion density gradient r]i (although we do not investigate in great detail the effects 
of changing rji). We find that with rjs = 0, A^max ~ 1 when g/e ~ 7 and that A^max in- 
creases monotonically with q/e. We find that as rjs — )■ 0, A^max tends to a constant 
determined by q/e, and that as rjs becomes large, A^max increases monotonically 
with it. 

Regarding subcritical turbulence, a prediction can be made from the results of 
Chapter [131 If we assume that in order for subcritical turbulence to be sustained 
nonlinearly, fluctuations must grow somewhat before they decay in order to pro- 
vide energy to be scattered into still-growing fluctuations, and if we assume that 
therefore A^max ^ A^c for subcritical turbulence to be sustained, where Nc, the crit- 
ical amplification exponent, is a number of order unity, the results predict that 
the PVG instability alone cannot drive turbulence for q/e below a certain value 
determined by Nc (for example, if Nc ~ 1, the PVG drive would be ineffective for 



In Chapter [141, inspired by this prediction, we use full nonlinear simulations 
to map out, as a function of the three key parameters, 7^;, q/e and k, the surface 
within that parameter space which divides the regions where subcritical turbu- 
lence can and cannot be sustained: the zero-turbulence manifold. As in Part [Till 
we take the Cyclone Base Case parameter regime jjzsl], i.e., concentric circular 
flux surfaces with e = 0.18, inverse ion density scale length R/Ln = 2.2 and ion 

^ The time varying radial wavenumber is defined in equation | |6.4|| . 



decays. 




(12.2) 



g/e < 7). 
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to electron temperature ratio Tj/Tg = 1 and magnetic shear s = 00 The ratio q/e 
is varied by varying q alone. The resolution of all simulations in Chapter [T4l was 
128 X 128 X 40 X 28 X 8 (poloidal, radial, parallel, pitch angle, energy). Note that 
relatively high parallel resolution was used to resolve the PVG modes, for the 
reasons given above. 

We cannot conclusively determine that there is a value oi q/e below which the 
PVG is unable to drive subcritical turbulence, but we are able to conclude that 
for q/e < 7, it is stable for all experimentally relevant values of the toroidal flow 
shear. We find that, as a consequence of this, at low q/e the maximum critical 
temperature can be very high: exceeding 20 for values oi q/e < 6, that is, com- 
parable to or hi gher than those experimentally observed for internal transport 



barriers 1152 



112(1 in advanced high power operating regimes. 



We also discover that the zero-turbulence manifold can be parameterised as 
Nma.y.ilEi Q'/c, = where A^max IS calculatcd from linear theory, but that 
is a function of 7^;. We are unable to explain this dependence. Nonetheless, the 
parameterisation would allow the full zero-turbulence manifold to be calculated 
using linear simulations once Nci^s) had been determined using a set of nonlin- 
ear simulations at a single value oi q/e. reducing the number of parameters that 
would be needed in the nonlinear calculation by one. 



* A temperature ratio of 1 is appropriate for both lower power ar i d fu ture reactor-like condi- 
tions, but not for high-performance shots in current devices Il09l lid. 111 1. 
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Chapter 13 
Subcritical Fluctuations 



This chapter is largely taken from Ref. nlOal . to which the author of this work 
was a major contributor. 

13.1 Introduction 

Where the magnetic shear is zero and the flow shear is finite, a fluctuation 
may only grow transiently before decaying. Nevertheless, it is possible, across 
a wide range of parameters, for this transient growth to be vigourous enough 
to allow subcritical turbulence to be sustained via nonlinear interactions. It is 
desirable, therefore, to be able to characterise the strength and duration of this 
transient growth, and to see whether it is possible to predict from that strength 
and duration whether or not turbulence could be sustained given a particular set 
of plasma parameters. 

13.2 The Flying Slab 

We will examine this transient growth using simplified slab geometry. This 
allows progress to be made analytically in certain limits while still including the 
main physical features of the ITG drive, the PVG drive and the perpendicular 
shear suppression. We will check and supplement these analytical results using 
the code Ast roGK , which is a faster, simplified version of GS2 designed for solv- 
ing problems in slab geometry. The gyrokinetic equation in slab geometry can be 
simply obtained from the full axisymmetric equation by taking circular geometry 
and letting i?— )■ oo, r — )■ oo, r/i?— )■ (which eliminates particle trapping) and 
Ruq — )• m, a constant. We further set as = 0. 



13.2. The Flying Slab 



We are, in effect, looking at a finite piece of an infinite torus; within that finite 
piece, all field lines are locally flat. The only geometric parameter which remains 
is s, which we set to 0. We note that the gyrokinetic equation (|3.45|) is in the frame 
rotating with velocity uq', our slab is therefore travelling at a finite speed u around 
the infinite torus, and so we refer to our model as the flying slab approximation. 



13.2.1 The Gyrokinetic Equation 



Starting from equation (|3.45|) , taking the limits given above, and defining 

d_ , d 

dz 



and 



B 

■^slab — j~) -^y 



^_Be dRu B^ 
^ = -B^ - 5"^^' 



(13.1) 
(13.2) 



(13.3) 



we obtain the flying slab gyrokinetic equation (which we will solve for ions only): 



7T7 + ^^slab'S'T^ ) hi 

at ay ^ 



+ 



(I+^^^^^^I) ^^^^^ 

2^11 ? 



Lfii 



££ _ 3 
Ti 2 



1 

Lti 



PiVth 



Rs 



dy 



(13.4) 



where = ZiCip/Ti is the normalised potentially] Once again we take the electro- 
static limit with a Boltzmann electron response, so that the system is closed by 
the quasi-neutrality condition: 



1 + 



Hi 



(fv{h) 



r' 



(13.5) 



where r = Ti/T^. For the remainder of this chapter we suppress the species 
indices and the slab label. 



13.2.2 Case of non-zero magnetic shear 

As a digression, we note that including a (locally) constant linear magnetic 
shear into the problem amounts to replacing in (|13.4[) 

„ dh ( ^ v\\\ dh 

Sx—^ [S+-^)x—, (13.6) 
dy V LsJ dy 



^ It should be noted that in this chapter we are dealing with a much simpler system than 
elsewhere in this work and hence choose a different system of normalisations. 
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13.2. The Flying Slab 



where Lg is the scale length associated with the magnetic shear. This appears 
to introduce complications as we now have an "effective shear" that depends 
on the particle velocity v\\. However, since the size of v\\ is constrained by 
the Maxwellian equilibrium distribution, this term can be neglected provided 
Ms = SLg/vth^ ^ 1. Under this assumption, the theory developed below applies 
without modification. We note that the "shear Mach number" Mg is precisely the 
parameter that is known to control linear stability and transient amplification of 
in the ITG-PVG-driven plasmas in the fluid (collisional) limit |49fl. 



13.2.3 Shearing frame 

The next step — standard in treatments of systems with linear shear — is 
to make a variable transformation (t, r) — )• (t', r') that removes the shear terms 

(Sxd/dy): 

t' = t, x' = X, y' = y — Sxt, z = z, (13.7) 

and similarly for (t, R) — {t', R'). The Fourier transform can then be performed 
in the primed variables, so 

k' k' 

where = k'^ — Skyt', ky = k'y, and k\\ = k'^i^ (we denote k\\ = k^,). As usual 
in the gyrokinetic theory, working in Fourier space allows us to compute the 
gyroaverages in terms of Bessel functions: 

(^)r = U<t'))^At'y'''-'^\ {h)r = Yl U<t'))hAty'-'\ (13.9) 



where a(t') = k^v^/n^ = [v^/ni]^ {k',- Sk'yt'f + k'^ = {k'yV^/ni)^/l + SH"\ 
and we have shifted the origin of time: t" = t' — S'^k'^/k'y. 

Finally, we rewrite the gyrokinetic system (|13.4|) - (|13.5|) in the new variables 
(t", k'). Since we are interested only in the linear problem here, we will drop the 
nonlinearity. We also suppress all primes in the variables. The result is 

-^ + ik\\v\\K = 

{-^-i uj*+(^ (^*Vi - —- Skypi \ FoJo{a{t))(^k, (13.10) 

(l + ^) = ^ y d\ Ma{t))h^, (13.11) 

where = kyPiVthj2Ln is the drift frequency and rji = Ln/ L^; the argument of 
the Bessel function is a{t) = {v±/vth JkyPi\^l + SH'^. 
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13.3. Solution for the case of strong shear 



13.2.4 Integral equation for the linearised problem 



We integrate (|13.10|) with respect to time, assume the initial fluctuation ampli- 
tude small compared to values to which it will grow during the subsequent time 
evolution, rescale time m \vthj — > t, denote At = t — t', use (|13.11|) , in which the 
velocity integrals involving the Maxwellian Fq = n e~^''^^^'^= / (vrf^j 
the usual way, and obtain finally the following integral equation for (fi{t): 

ft 



I (vrf^\ )^/^ are done in 



9 - 



l + r^, A(A,A')-1- 



4 



ro(A,A')^(t-At), 



(13.12) 



where = iij^l\k\\\vw,^ = kyPi/2\k\\\Ln is the normalised drift frequency, ujs 
SkyPi/k\\VtYi, the normalised shear parameter, and 

A + A' , /rryMy^^ 



ro(A, A') = e-(^+^')/^/o(v^), A(A, A') 



^AA 



^AA' 



(13.13) 



where X{t) = {k^p^ + uln/2, X' = A(t') = (/t^pf + ult'^)/2, and Jq and h are 



modified Bessel functions of the first kind. 



Equation (|13.12|) is the master equation for the linear time evolution of the 
plasma fluctuations driven by the ITG (the u)^7]i term) and the PVG (the {q/e)ujs 
term). 



13.3 Solution for the case of strong shear 

We will first consider the maximally simplified case of pure PVG drive (strong 
shear). This is a good quantitative approximation to the general case if uo^,r]iUo^ <ti 
(g/e)a}5, which in terms of the basic dimensional parameters of the problem trans- 
lates into 

^ » ^, ^. (13.14) 

This is the regime into which the plasma is pushed as the flow shear is increased 
— under certain conditions, the transition can occur abruptly, via a transport 
bifurcation (see PartHIIl). Besides being, therefore, physically the most interesting, 
this limit also has the advantage of particular analytical transparency (the more 
general case including ITG will be considered in Sec. 113.4)) . 

Thus, neglecting all terms that contain and r]i, (|13.12|) becomes 

(l + ^ - ro(A, A)) m = 2 y dAt At e-^'"l' c:;^ - l) ro(A, A')^(t - At). 

° (13.15) 
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13.3.1 Short-time limit: the PVG instability 

Let us first consider the case in which the velocity shear is unimportant except 
for the PVG drive, i.e., we can approximate X ^ X' ^ and so there is no 

time dependence in the Bessel functions in (|13.15|) . We would also like to be able 
to assume t ^ 1 so that the time integration in (|13.15|) can be extended to oo. 
Formally this limiting case is achieved by ordering kyPi ~ 1 and 1 <^ t <^ 1 /cus 
(in dimensional terms, this is equivalent to St <^ 1, k\\vthj ^ 1 and k\\vthj S ^ 
1). Physically, this regime will be realised in the initial stage of evolution of the 
fluctuations and, as we are about to see, will require a large enough q/e. 

Under this ordering, we can seek solutions to (|13.15|) in the form Lp{t) = 
(po exp(— iwt), where u) = uj/\k^\ |fths is the nondimensionalised complex frequency 
and u! its dimensional counterpart. The time integral in (|13.15D can be expressed 
in terms of the plasma dispersion function, which satisfies 1113(1 



/■oo 

Z{u) = i / dAte'^^'-^''/\ Z\uj) = -2[1 + uZ{u)]. 
Jo 



(13.16) 



With the aid of these formulae, (|13.12[) is readily converted into a transcendental 
equation for u: 

1 + t/Z ^ (q_ 



(-ws-l) [l + wZ(u;)], (13.17) 



ro(A) Ve 
where ro(A) = e-^/o(A), A = A;^pf/2 and Qs = Skyp,/k\\Vti,,. 

Equation (|13.17|) is simply the dispersion relation for the ion acoustic wave 
modified by the PVG drive term. This point is probably best illustrated by 
considering the cold-ion/long-wavelength limit r ^ 1, u; = a;/|A;|| |t>ths ^ 1/ 
A = klp^/2 < 1. Then ro(A) ^ 1, 1 + ojZ{uj) ^ -1/2lo^ + iu^ e-^\ and so, 
restoring dimensions in (|13.17|) , we get 

c^-^fl-^a;,) ^ a;-±fc,|cJl-^A^V'', (13.18) 



2r V e / " V e k\\Vt\,, 

where in the last expression, we have restored dimensions and denoted = 
(Z/2r)^/^fth^ = {ZTe/niiY/'^, the sound speed. When {q/e)us is sufficiently large, 
the sound wave is destabilised and turns into the PVG instability. Note that it 
loses its real frequency in this transition. 

Some further (elementary) ana lytical considerations of the PVG instability are 
given in Appendix B.l of Ref . |108I1 . Here it will suffice to notice that if the (dimen- 
sional) growth rate 7 = Imw is scaled by qS/e and k\\ by iqS/e)/vth,, their mutual 
dependence is universal for all values of the velocity shear or g/e, namely, 

7=^/fv..^V (13.19) 
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Figure 13.1: The PVG Instability 



Left panel: Growth rate 7 = Imcj (normalised by qS/t) vs. k\\Vths/iqS/e). Here 
t/Z = 1, and the two curves are for kyPi = 0.1 (red) and kyPi = 1 (brown). 
The growth rate becomes independent of ky for kyPi ^ 1 (see end of Appendix 
B.l of Ref. jlOSIl ) and all curves for large kyPi are very close to each other and 
similar in shape to kyPi = 1 (see right panel). The mode has no real frequency 
for k\\Vths / {qS / e) < ^yP*/ foi" kpths / (qS / e) > kyPi, it turns into a damped sound 
wave: the corresponding frequencies and the damping rate are shown only for 
kyPi = 0.1, as thin blue (dashed) and red (solid) lines, respectively. The discrete 
points show growth rates calculated by direct linear numerical simulation us- 
ing the gyrokinetic code AstroGK [72]. Right panel: Contour plot of 7/(gS'/e) 
vs. kpths/ilS/e) and kyPi. Only positive values are plotted, black m eans 7 < 0. 
The red curve shows the stability boundary (equation (B.l) of Ref. 1108(1 ). The 
white line shows the boundary (equation (B.2) of Ref. [108]) between PVG modes 
(above it) and sound waves (below it). 



Once this rescaling is done, the dispersion relation (|13.17|) no longer contains 
any parameters (except r/Z, which we can safely take to be order unity). The 
growth rate, obtained via numerical solution of (|13.17[) with r/Z = 1, is plotted 
in Fig. 113.11 The maximum growth rate is 7max ~ 0.10(g5'/e). This peak value is 
reached when kyPi ^ 1.0 and k\\ ^ 0.10(gS'/e)/fth^, i.e., at {q/t)uos ~ 10. 

The conclusion is that, at least in the initial stage of their evolution, plasma 
fluctuations in a significant part of th e wavenumber space {ky and k\\ are con- 
strained by (equation (B.l) of Ref. |108|] )) are amplified by the PVG. This amplifi- 
cation does not, however, go on for a long time as the approximation we adopted 
to derive the dispersion relation (|13.17[) breaks down when ust ~ kyPi (or St ~ 1 
if the dimensional units of time are restored) — this gives 7t ~ 0.1 (g/e), so, real- 
istically, after barely one exponentiation. The key question is what happens after 
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that. We will see shortly that all modes will eventually decay and that the fastest 
initially growing modes are in fact not quite the ones that will grow the longest 
or get maximally amplified. 



13.3.2 Long-time limit: transient growth 

Let us now investigate the long-time limit, cost ^ 1, uJst ^ kyPi (or, in di- 
mensional form, St ^ 1, kyPiSt ^ 1). In this limit, the kernel involving the 
Bessel function in (|13.15|) simplifies considerably: we have A ^ 1, 

A' f« - At) 2/2 > 1 and so 



ro(A,A') 



-(^A-^/V)V2 



-(D|AtV4 



(13.20) 



Working to the lowest nontrivial order in 1/t, we can now rewrite (|13.15|) as fol- 
lows 



1 + 



r 



\ujs\tLp{t) 



1 







^y,-orrv^/ 20F 

We will seek a solution to this equation in the form 



l)^(t-At). (13.21) 



(13.22) 



where ^{t) = 7(t)/| A;|| \vth^ is the effective time-dependent growth rate (nondimen- 
sionalised) and 'y{t) the dimensional version of it (remember that time is scaled 
by |A;| |fthj. Because of the exponential in the kernel, the memory of the time- 
history integral in the right-hand side of (|13.21|) is limited, so At ^ t and we will 
be able to make progress by expanding 

At^ 



if{t — At) = yjQ exp 



/ cit'7(t') - At7(t) + 
Jo 



.f(t) + ... 



(13.23) 



We will assume that this expansion can be truncated; the resulting solution will 
indeed turn out to satisfy 7'(t) ^ 1, with all higher-order terms even smaller. 
Substituting (|13.23|) into (|13.21|) , we obtain an implicit transcendental equation 
for 7(t): 

(l + l\\Q^\t = A^ /"dAt Ate-^*^W-(^i+i)AtV4 - iV (13.24) 

This equation can be written in a compact form by invoking once again the 
plasma dispersion function (|13.16|) : denoting 7(t) = 7(t) / ^Juj^ + l, we get 

(l + ^) v^(4 + l)\^s\t =[-^cos- l) [1 + t^mit^m (13.25) 

— effectively, a time-dependent dispersion relation, reminiscent of the PVG dis- 
persion relation (|13.17|) . 
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Transient growth. 

First of all, since the effective time-dependent growth rate 7(t) appears with 
a negative sign in the exponential in (|13.24[) , it is clear that as time increases, (the 
real part of) 7(t) must decrease and indeed go negative because the right-hand 
side has to keep up with the increasing left-hand side. Therefore, fluctuations 
will eventually decay. However, if Re7(t) is positive for some significant initial 
period of time, there can be a substantial transient amplification. 

We can determine the time to when the transient growth ends by setting 
Re 7(^0) = in (|13.25)) . We immediately find that Im7(to) = as well and that 



r/Z)y/^{Ql + l)\us\' 



(13.26) 



The dependence of to on ky and k\\ — via ojs = SkyPi/k^^Vths arid via the time nor- 
malisation factor of |fc||fths I — tells us which modes grow longest. The interesting 
question, however, is rather which modes get maximally amplified during this 
transient growth. 
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Figure 13.2: Time evolution of the effective growth rate: the red (bold) line is 
7(t) = l{t)/\k\\\vths\/^s + 1 obtained as a numerical solution of (|13.27|) and plot- 
ted vs. t/to (the time axis is logarithmic in base 10). The black (thin) lines are the 
asymptotics (|13.28|) and (|13.33|) . The discrete points show the time evolution of 
the effective growth rate obtained in a direct linear numerical simulation using 
the gyrokinetic code AstroGK [72] with = 1/Lt = 0, q/e = 50, r/Z = 1, 

kyPi = 1, k\ivthj S = 0.5. The short- time-limit PVG growth rate for this case 
(obtained by solving (|13.17[) ) is shown as a dotted horizontal line. The time is 



normalised using (|13.26|) for to- 
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Maximal amplification. 

The total amplification factor is given by , where = Jq" dt'j{t) is the num- 
ber of exponentiations experienced by the mode during its growth period. In 
order to determine this, we need to know the time evolution of ^{t) up to t = to- 
Using (|13.26)) , it is convenient to rewrite (|13.25|) as follows 



- = l + i^Z{i^). (13.27) 

When t <^ to, the solution is found by expanding the plasma dispersion function 
in 7 > 1: 

to 272 ^ ^ V 2t 
This asymptotic is not valid when t approaches to. More generally, (|13.27[) has 
a solution 7 = 7(t/to)/ whose functional form is independent of any parameters 
of the problem. It is plotted in Fig. 113.21 together with the asymptotic (|13.28|) 



~. ' (13.28) 



and with a direct numerical solution showing how the transition between the 
short-time (Sec. 113.3. 1|) and long-time limits occurs. The amplification exponent 



is easily found: 

N = r dm = to^^i f d^m - 0.45 7r-^j^^%r=T7- 

Jo ^ Jo {1 + T/Z)y/C0g + l\UJs 



where we have used (|13.26|) for to arid computed the integral numerically. It is 



clear from (|13.28|) that the integral converges on its lower limit and is not domi- 



nated by it, so it does not matter that we cannot technically use (|13.27|) for short 



times (as confirmed by comparison with a direct numerical calculation in the 
right panel of Fig. 113. 3|) . 

According to (|13.29)) , N depends on both wavenumbers via cus only (a plot of 
this dependence will be given in Fig. 113. 6|| . Assuming q/e ^ 1, we find that the 
amplification exponent is maximised for cus ~ {^/qY^^, giving 

iVmax « 0.45 for kyp,^(^^y (13.30) 

These results are illustrated in the left panel of Fig. 113.31 The amplification time 
for the maximally amplified modes identified in (|13.30|) is, from (|13.26)l , 

, ^, q/e jq/er/' 



^Note that this is well outside of the wavenumber domain populated by the damped sound 
waves, LOS < e/q (see Appendix B.l of Ref. [108.1 ). 
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Figure 13.3: The PVG Amplification Exponent 

heft panel: The amplification exponent N vs. k± and /cy in the same AstroGK 
simulation as used to produce the discrete points in Fig. 113.21 The straight line 
shows the relationship between the wavenumbers given by the second formula 
in (|13.30|) . Right panel: The maximum amplification exponent A^max [N given by 
(|13.29[) , maximised with respect to us] vs. q/e. The dotted lines show the q/e:^ 1 
asymptotic [see (|13.30|) 1 and the g/e ^ 1 asymptotic, the latter straightforwardly 
obtained from (IT3^ : us ~ 2(e/g), iV„,ax ~ 0.11(g/e)V(l + r/Z) (but this is purely 
formal because the long-time conditions St^ 1 and kyPiSto ^ 1 will be bro- 
ken in this regime, except for extremely long wavelengths). The discrete points 
show A^max obtained via an AstroGK numerical parameter scan: k\\Vths/ S = 0-5/ 
varying k±pi and holding all other parameters fixed as in Fig. 113.21 (note that the 
asymptotic results do not in fact depend on k\\ or S, although the quality of the 
long-time asymptotic does). 



where we have restored dimensions to make explicit the dependence of to ori /cy . 
Thus, we have learned that an entire family of modes, characterised by a par- 



ticular (linear) relationship between ky and A;||, given in (|13.30|) , will eventually 
enjoy the same net amplification, even though, as follows from the results of Sec. 
113. 3. 1[ they were not the fastest initially growing modes and some within this 
equally amplified family started off growing more slowly than others or even de- 
caying. The more slowly growing modes are the longer-wavelength ones and, 
according to (|13.31|) , they compensate for their sluggishness with longer growth 
times. 

Note that (|13.31|) confirms that the long-time limit is analytically reasonable 
because for large q/e, (|13.31|) formally satisfies kyPiSto ~ {q/^Y^^ ^ 1 and also 
SIq ^ 1, provided A;||fth^/S' ^ q/e. The latter condition is marginally broken by 
the fastest initially growing modes: indeed, in Sec. 113.3. 1[ we saw that they had 
^\Vths/ S ~ O.lg/e, so for them, St^ ~ 3, not really a large number. For all longer- 
wavelength modes, to is safely within the domain of validity of the long-time 
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limit. 

Note also that, as is shown at the end of Appendix B.l of Ref. |108l] , the time- 
dependent dispersion relation (|13.27|) and its consequences derived above can be 
obtained directly from the PVG-instability dispersion relation (|13.17|) simply by 



restoring its dependence, setting = Skyt, taking the short-wavelength and 
long-time limit (k^pi 3> 1, St ^ 1), and assuming e/q <^ us ^ I. This calculation 
underscores the fundamental simplicity of the physics of the transient amplifica- 
tion: perturbations initially destabilised by the PVG are eventually swept by the 
perpendicular velocity shear into a stable region of the wavenumber space. 

Limits on short and long wavelengths. 

We have seen that modes with parallel wavenumbers up to /cyt'ths /S" ~ g/e can 
be transiently amplified. From (|13.30[) , we conclude that of these, the maximally 
amplified ones will have perpendicular wavenumbers up to kyPi < (g/e)^/^, i.e., 
kyPi can be relatively large — unlike in the short-time limit treated in Sec. 113.3. 1[ 
where the modes with kyPi ~ 1 grew the fastest (although large kyPi were also 
unstable). 

It should be understood that, while there is no ultraviolet cutoff in our theory 
that would limit the wavenumbers of the growing modes (in either direction), 
such a cutoff does of course exist in any real system. In the parallel direction, 
those fc|| that were strongly damped in the short- time limit (see equation (B.l) of 
Ref. ilOSO) are unlikely to recover in the long-time limit. In the perpendicular 
direction, the cutoff in kyPi will come from the collisional damping, which, in 
gyrokinetics, contains a spatial diffusion (see, e.g., izd]), and from the electron 
Landau damping, which we have lost by using the Boltzmann electron response 
(see (|13.5[) ) and which should wipe out large ky and k^\. 



On the infrared (long-wavelength) side of the spectrum, we have no cutoffs 
either. In a slab, these would be provided by the dimensions of the periodic 
box. In a real plasma, the cutoffs are set by the scales at which the system can 
no longer be considered homogeneous (in a tokamak, these are the equilibrium- 
gradient scale lengths and the minor radius for the perpendicular scales and the 
connection length qR for the parallel scales). 

Significant amplification threshold. 

If we maximise (|13.29|) without assuming q/e ^ 1 (with the caveat that the 
long-time limit asymptotics are at best marginally valid then), we obtain a more 
general curve than (|13.30)) , plotted in the right panel of Fig. 113.31 We may define 
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a critical threshold for significant amplification: A'n^ax = 1 when g/e ^ 7. The role 
of this threshold will be discussed in Chapter [T4l 

Long-time decay. 

Finally, we obtain the long-time asymptotic decay law. Let us seek a solution 
of (|13.27|) such that t > to and 7 < -1. Then 



^^2v/^|7|e^' 7^-^ln^. (13.32) 

Since 7 is only root-logarithmically large, the quality of this asymptotic is rather 
poor. If we insist on a more precise decay law, we can retain small corrections in 
(|13.32[) and get what turns out, upon a numerical test, to be a reasonably good 
approximation (see Fig. I13.2|) : 




^/(2v^^o) (13.33) 



For the maximally amplified modes (see (|13.30|) ), the dimensional damping rate 
is 7(t) ^ \k\\\vth^^{t) with to given by (|13.31|) . This tells us is that the decay is just 



slightly faster than exponential at the rate of order S. The longest-wavelength 
modes decay the slowest, after having being amplified the longest. 

13.4 Solution including ITG 

Let us now generalise the results obtained in Sec. ll3.3l to include non-zero (i.e., 
non-negligible) density and temperature gradients. This means that we restore 



the terms involving oj^ and r/j in the general integral equation (|13.12[) . 



13.4.1 Short-time limit: the ITG-PVG dispersion relation 

The short- time limit, introduced at the beginning of Sec. 113.3.11 for the case 
of pure PVG, is treated in an analogous fashion for the general ITG-PVG case. 
An analysis of the solutions of the resulting dispersion relation is useful in that 
its results assist physical intuition in ways relevant for some of the forthcoming 
discussion, but it is not strictly necessary for us to have them in order to work out 
how transient growth happens in the presence of the ITG drive. For this analysis 



therefore we refer the interested reader to Appendix B.2 of Ref. |jl08(]. 
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13.4.2 Long-time limit 



We now continue in the same vein as in Sec. |13.3!2l and consider the long-time 
limit {St ^ 1, kyPiSt ^ 1), in which we can simplify the kernels involving the 
Bessel functions in (|13.12)) by using (|13.20|) and also 



A(A,AO 



1 (v^-v^)2 1 ul/\t^ 



2 2 2 4 

This allows us to rewrite (|13.12[) in the form that generalises (|13.21|) : 



(13.34) 



1 + ^) \i^smt) = ^ 







°°rfA«c-<"l+'>'^'"''«|(ii.s-l)^ 

1_J1^M±1)^)1U„_A«). (13.35) 



As in Sec. 113.3. 2[ we seek a solutions to this equation in the form (|13.22|l , taking 
At <^ t and expanding the delayed potential under the integral according to 
(|13.23|) . The result is the generalised form of (|13.24|) : 



r 



1 + \^s\t 



dAt e-At7{t)-{4+i)AtV4 ^ 



At 



lOJ* 



. (13.36) 



Using again the plasma dispersion function (|13.16[) to express the time integrals 
and introducing the complex scaled frequency a)(t) = 27 (t) / a/wI + 1, we get 



r 



y^7r(a;| + l)\ujs\t 



{q/e)us - 1 



- r]iUJ^uj{t) 



+ {r]i-l)u;,Z{u{t)), 



[l + u{t)Z{u{t))] 

(13.37) 



a time-dependent dispersion relation, which is the generalisation of (|13.25|) . 



Transient growth. 

The general argument that the real part of 7(t) (i.e., the effective time- 
dependent growth rate) must eventually decrease and so fluctuations will, in 
the end, decay, applies to (|13.36|) similarly to the way it did to (|13.24|) (see Sec. 
113.3. 2[) , although this decay need not (and, as we will see, will not) be monotonic. 
The time to when the transient growth ends is now determined as follows. Let 
Ima)(to) = and d)(to) = cDq (real!). Then from (|13.37[) taken at t = to, we find the 
real frequency coq by demanding that the imaginary part of the right-hand side 
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vanish — this means that the coefficient in front of Z{ujq) must be zero, because, 
for real Uq, ImZ(d)o) = ^¥^"'^0. This condition gives 

1 



- ~2 {(l/^)^S 



■ Uo - {rji - l)u^ = 0, 



(13.38) 



whence 



(13.39) 



Substituting this solution into the real part of (|13.37|) and taking advantage of the 
already enforced vanishing of the coefficient in front of Z(uq), we get 



{q/e)Os - 1 + ^/[{q/e)u;s - I? + - V^.) ^5 H + 1) 



(13.40) 



2{l + r/Z)^{ul + l)\us\ 

where we have replaced rjfu'^ = rjgUg/A with rjs = v^^JLtS a new parameter 
that measures the strength of the ITG drive relative to the velocity shear. Note 
that we picked the "— " mode in (|13.39|) because the "+" mode is not amplified 
(to < 0, assuming r^j > 1). Equation (|13.40[) is the generalisation of (|13.26[) , to 
which it manifestly reduces when rjs = Q and with which it shares the property 
that the transient growth time depends on ky and only via ojs arid the time 
normalisation factor |/i;ii|fth • 



General dispersion relation. 

We can now recast the general ITG-PVG case in a form that shows explicitly 
how it reduces to the case of pure PVG drive studied in Sec. 113.3.21 First we note 
that the transient growth termination time (|13.40|) can be rewritten as 



.(PVG) 



X 



{q/e)ujs-l 



(13.41) 



where tg^^^-* is given by (|13.26|) , rjsUJs = 2riiUJ^, rjs = Viy^JLtS and a 



sgn[(g/e)a;5' — 1]!^ Then the time-dependent dispersion relation (|13.37[) can be 
manipulated into the following form: 



1 + (tJ1 + 



1 

1 — 



f = (2 - X cD) [1 + Z{Co)] +(l--]^Z{Co). (13.42) 



The analogous equation for the case of pure PVG drive, (|13.27[) , is recovered when 
X <S 1, which means 775 ^ q/e (cf. (|13.14|) ) and "qsuis ^1- In this limit, the 

^If (g/e)ws < 1, 4^^^^ < and X < 0, but the ITG mode can still have transient growth, 
to > 0. 
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behaviour of fluctuations in the presence of both PVG and ITG drives is well 
described by the results of Sec. 113.3.21 Before discussing the general case, it is 
useful to consider the opposite extreme of weak velocity shear. 

13.4.3 Case of weak shear 

Let rjs q/e and risUs ^ 1, so x ^ 1 (note that the same limit is also achieved 
for ids — )■ oo). Then the x dependence falls out of (|13.42)) : 



- = -u[l + CdZ{Co)]+ Z{Co), (13.43) 

^0 

where we have discarded the I/t^j terms by assuming r]i ^ 1. The transient 
growth termination time in this limit is, from (|13.41|) , 



to = -X tr""^ = (13.44) 

When t <^ to, we find the solution of (|13.43|) by expanding in a) ^ 1. It turns out 
that it consists of a large real frequency and an exponentially small growth rate: 
(|13.43D becomes 

l.__L_.V^..,-- ^ ,._| + .2V^(|)'e-<'"^"'. (13.45) 

More generally, for finite values oit/to, the solution of (|13.43|) is d) = cI;(t/to), with 
a functional form independent of the parameters of the problem. This solution 
it plotted in Fig. 113.41 It turns out that at t ^ 0.15 to/ the growth rate increases 
sharply, reaches a finite maximum and then decreases towards zero, which it 
reaches at t = to, whereupon growth turns to decayE 

Thus, there is a period of strong transient amplification, which lasts for a finite 
fraction of time to. The amplification exponent is 

N = jJdtj{t)=to^ul + lj^ d^lmuiO- 0.057^^^, (13.46) 

where we have used (|13.44[) and calculated the value of the integral under the 
curve in Fig. 113.41 numerically (note that since the growth rate is exponentially 
small at t <^ to, the precise lower integration limit is irrelevant). Remarkably, 

■^This implies that perturbations first grow due to the ITG-PVG instability (at St ^ 1), then 
slow down to exponentially small growth rates, then (at St ^ 1) grow vigorously again before 
finally starting to decay at t = to. The intermediate period of virtually zero growth, which is 
a feature both of the weak-shear regime and of the general case (see Sec. 113.4.4b and may ap- 
pear strange at first glance, can be traced to the dependence of ITG -PVG growth rates at long 
parallel wavelengths — this is explained in Appendix B.2 of Ref. jlOSll . 
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Figure 13.4: Time evolution of the effective growth rate and frequency: the red 



(upper bold) line is 7(t) 
is Reu{t)/\k 



Irnw(t)/|A;|||t'th^A/w|~TT and the blue (lower bold) line 



^yuJg + 1, both obtained as a numerical solution of (|13.43[) and 



Vths V 

plotted vs. t/to (the time axis is logarithmic lo). The black (thin) lines show the 
gro wth r ate and frequency given by the asymptotics (|13.45|) and equation (C.9) of 
Ref. jlOSl] (the latter taken in the limit x — )• oo and r]i oo). 



unlike in the case of the PVG drive (see (|13.30|) ), the amplification exponent has 



no wavenumber dependence at all. Also unlike in the PVG case, it does depend 
on the shear and on the temperature gradient: N (xrjs = Vths/LxS. 

To recapitulate, we have found that, at low velocity shear, all modes are ampli- 
fied by a large (and the same) factor before decaying eventually. Their transient 
amplification time is given by (|13.44|) . Restoring dimensions, (|13.44|) and (|13.46|) 



are 



to- W(£t£)^^^ at ^ 0.057 (13.47) 

2{1 + T/Z)Jn{S^klpj + klvlJ 1 + ^/^ 



The transient growth lasts for a very long time at low S and the longest-growing 
modes are the long-wavelength ones. The limit S* — is singular in the sense that 
for arbitrarily small but non-zero S all modes eventually decay, while for S = 0, 
the indefinitely growing linear ITG instability is recovered (to = oo, N = oo). 
We have already made the point (in Sec. ll3.3.2|) that while our theory does not 



limit the transiently growing wavenumbers from above, a fuller description of 
the plasma will. 
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t 

to 



Figure 13.5: Effective normalised growth rates lraLd{t/to) (red, top) and fre- 
quencies Re a) (t /to) (blue, bottom) obtained via numerical solution of (|13.42[) with 
rji = 5 and Y = 0.1, 1, 2, 10 (from top/lighter to bottom /darker curves). See Fig. 
CI of Ref . jlOSl] for a more detailed depiction of the x = 1 case. 



13.4.4 Case of finite shear 

In the intermediate regimes between large and small x (i-e-/ weak and strong 
shear), the solutions of (|13.42|) transit from the weak-shear form described in Sec. 
113.4.31 to the pure-PVG case treated in Sec. 113.3.21 Figure 113.51 shows the time- 
dependent growth rates and frequencies for several values of x- As x decreases 
(i.e., S increases), the peak of the growth rate moves further into the past and 
the growth rate asymptotes to the pure-PVG case (Fig. |13.2|) . It is not hard to 
convince oneself analytically that this is indeed what ought to happen. F urthe r 



asymptotic considerations on this subject are given in Appendix C of Ref. |108|1 . 
The long-time decay asymptotic is also derived there (it is exactly analogous to 
that found in Sec. I13.3.2D . 

Similarly to our previous calculations, the amplification exponent is 

N{us) = to{Qs)\Jujl + l dilmQ{i,x{oJs)). (13.48) 

where Co is the solution of (|13.42|) . The wavenumber dependence enters via the 
ujs dependence of to arid of x (see (|13.41[) ). The numerically computed amplifica- 



tion exponent as a function of uos and of r]s is plotted in the left panel of Fig. 113.61 
for rji = 5 and q/e = 10 (these are representative of the values encountered in 



the more realistic numerical studies of tokamak transport |45l,l46l,l88|1 ). The pure- 
PVG case treated in Sec. 113.3.21 remains a good approximation up to values of rjs 
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Figure 13.6: The PVG-ITG Amplification Exponent 



The amplification exponent N, given by (|13.48|) with u the solution of (|13.42|) for 
r^. = 5, q/e = 10 and r/Z = 1. Left panel: N vs. ^5 for r]s = 1, 2, 5, 10, 20, 50, 100, 200 
(from bottom to top curves, in darkening shades of red). The numerical results 
for the same cases are shown in Fig. 113. 7[ The r]s = case (pure PVG drive) 
was virtually indistinguishable from r]s = I when plotted (not shown here). The 
lowest-r^s behaviour is well described by (|13.29|) , the highest-r^s' by (|13.46|) (in 
the latter case, except for corrections associated with finite rji, which are easy to 
compute if they are required). Right panel: The maximal amplification exponent 
N vs. the normalised shear (77^^ = SLt /vths)- The maximum N(tus) is reached 
at Us ~ 0.54, independently of rjs. The dotted lines are the asymptotics (|13.30|) 
and (|13.46|) . The finite offsets between the asymptotics and the exact curve are 
due to the fact that the asymptotics were calculated in the limits r/j ^ 1 and 
g/e ^ 1, while for the exact solution we used relatively moderate values of these 
parameters. The discrete points show A^max obtained via an AstroGK numerical 
parameter scan varying k±pi and 5* (i.e., r]s) while holding k\\LT = 0.02 and the 
rest of the parameters fixed at the same values as quoted above. 



of order 10. After that, there is a transition towards the weak-shear limit (Sec. 
I13.4.3[) , accompanied by the loss of wavenumber dependence as rjs is increased 
to values of order 100. The constant of proportionality between ky and k\\ for the 
maximally amplified modes (i.e., at which N is maximised) does not appear 
to depend on r]s, although at large rjs, the maximum is increasingly weak. The 
maximal amplification exponent is plotted in the right panel of Fig. 113.61 It is per- 
haps worth pointing out the qualitative similarity between this plot and figure 1 
of Ref . ||46^, obtained from gyrokinetic simulations in full tokamak geometry. 

Finally, the amplification exponent as a function of ky and /cj , obtained in di- 
rect (linear) numerical simulations, is shown in Fig. 113.71 (the parameters are the 
same as in the "theoretical" Fig. 113. 6[ left panel). The transition from the PVG 
curve (|13.30[) to the flat wavenumber dependence (|13.46[) is manifest, as are the 
limits of applicability of our approximations in the wavenumber space. Note 
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13.5. Qualitative summary of the linear results 



the different normalisation of the parallel wavenumber here {k\\LT, characteristic 
of ITG) compared to the left panel of Fig. 113.31 (A;||ftv„ /5', characteristic of PVG). 
Hence the drift towards higher k\\LT as r]s = Vthi/SLx decreases towards the 
PVG-dominated regime, where the parallel scale of maximally amplified modes 
is set by the shear rather than the temperature gradient. 



13.5 Qualitative summary of the linear results 

In a gyrokinetic plasma with radial gradients of temperature and parallel ve- 
locity, both gradients are sources of free energy and so will drive the growth 
of fluctuations (ITG and PVG instabilities). The typical growth rate is of order 
7 ~ Vths /Lt for ITG and 7 ~ qS/e for PVG (see (|13.19|) ), or the mean square of 
the two if they are comparable (see equation (B.IO) of Ref. [|108(1 ). Because the 
mean plasma velocity is toroidal, it always has both a parallel and a perpendic- 
ular component (the latter a factor of g/e smaller than the former). The shear in 
the perpendicular (E x B) velocity is stabilising and causes all modes to decay 
eventually, so the fluctuation growth is transient — it is always transient in the 
limit, considered here, of zero magnetic shear and it is transient for large enough 



velocity shear S when the magnetic shear is finite |49l. 145(1 . If the linear physics 
provides sufficiently vigorous and lasting amplification of finite initial perturba- 
tions, the system is able to sustain nonlinearly a saturated (subcritical) turbulent 
state (see Chapter [12)). Therefore, the interesting question is how much transient 
amplification should be expected to occur and on what time scale. 

In the preceding sections, we have addressed this question mathematically, 
with the results summarised by figures [T33lll3.6l and ll3.7l (see also (|13.30|) , (|13.31|) 
and (|13.47|) ). Very roughly, these results can be explained as follows. The effect 
of the perpendicular shear is to produce a secular increase with time of the ra- 
dial wavenumber, (t) ~ Skyt. When this becomes large enough, the instabil- 
ity is kill ed by Landau damping (see discussion at the end of Appendix B.l of 
Ref. ilOSD). If we estimate that this happens after to ~ S^^ (i.e., for kx{tQ)pi ~ 1, 
assuming kyPi ~ 1), we may conclude that initial perturbations will be amplified 
by a factor of , where the amplification exponent is 

AT ~ ~ for ITG and ~ - for PVG. (13.49) 

Thus, the shear quenches the ITG amplification — but N cannot fall below the 
shear-independent level associated with the PVG (Fig. 113. 6[ right panel). This 
is indeed the case (see (|13.30|) and (|13.47|) ), although, strictly speaking, one has 
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to take into account the dependence of the quenching effect on the perpendic- 
ular and parallel wavenumbers — long-wavelength modes grow more slowly, 
but for a longer time; in the case of PVG, there is also a preferred relationship 
kypi ~ {e/qY/'^k\\Vtu/ S for the most strongly amplified modes (see Sec. 113.3.21 and 
Fig. 113. 7|) . While these wavenumber dependences are likely to be important in 
the analysis of the resulting turbulent state and the associated transport, they 
effectively cancel out in the expression for the amplification exponent (because 
7 oc kyPi, to oc l/kyPi at long wavelengths) and the results of the qualitative argu- 
ment that we have given hold true. 

It is instructive to compare these results with the conclusions of a long- 
wavelength fluid ITG-PVG theory presented in |49|1 (for the case of finite mag- 
netic shear). In that regime, perpendicular shear, by effectively increasing kx{t), 
also caused eventual damping of the fluctuations, but this time via coUisional 
viscosity. Therefore, to estimate the transient growth time to, one must set 
7 ~ iyiikl{to)pf ~ UiiS'^kypftl, where uu is the ion collision rate. Then, ignoring 
wavenumber dependences again, to oc 7^/^S'~^, so the amplification exponent is 

^3/2 ^ 

~ 7^0 OC ^ OC - for ITG and oc for PVG. (13.50) 

o o 

Thus, the E x B velocity shear again quenches the ITG instability, but once S 
is large enough for the PVG drive to take over, the amplification exponent actu- 
ally grows as v^, the result obtained rigorously by |49|1 — in contrast with the 
shear-independent ~ g/e that we have found in the kinetic regime. The prac- 
tical conclusion from this is that it should be easier to obtain states of reduced 
transport |45Ll46l,lll4|l in weakly coUisional, kinetic plasmas. 
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Chapter 14 
The Zero-Turbulence Manifold 



This chapter is published as Ref. /Ilia/ . 



At zero magnetic shear, the turbulence is subcritical for all nonzero values of 
the flow shear: there are no linearly unstable eigenmodes, and sustained turbu- 
lence is the result of nonlinear interaction between linear modes which grow only 
transiently before decaying. In Chapter [131, which studied this transient growth 
in slab geometry, it was demonstrated that at large velocity shears the maximal 
amplification exponent of a transiently growing perturbation before it decays is 
proportional to the ratio of the PVG to the perpendicular flow shear. In a torus, 
this quantity is equal to the ratio of the toroidal to poloidal magnetic field com- 
ponents, or q/e, where q is the magnetic safety factor and e is the inverse aspect 
ratio. Therefore, if we conjecture that a certain minimum amplification exponent 
is required for sustained turbulence, the results of Chapter [13] predict that there 
should be a value of g/e below which the PVG drive is rendered harmless. Below 
that value of q/e, it should be possible to maintain an arbitrarily high temper- 
ature gradient without triggering turbulent transport provided a high enough 
perpendicular flow shear can be achieved. 

In this chapter, motivated by the possibility of reduced transport at low values 
of q/t, we use nonlinear simulations to map out the zero-turbulence manifold, the 
surface in the parameter space that divides the regions where turbulent transport 
can and cannot be sustained. The parameter space we consider is {'Je, q/e, k). As 
ever, we set the magnetic shear to zero, the regi me we expe ct t o be most amenable 
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88 



114 [52, 



51 



load 



to turbulence quenching by shear flow |45l. 

We discover that reducing q/e is indeed uniformly beneficial to maintaining 
high temperature gradients in a turbulence-free regime, and that values of k can 
be achieved that are comparable to those experimentally observed for internal 



14.1. Finding the Boundary 



transport barriers | |52l,lll2l] . 

In the next sections, having presented our methodology, we will describe 
these results and discuss their physical underpinnings, as well as their impli- 
cations for confinement in a toroidal plasma. We will show that linear theory of 
subcritical fluctuations flOSl] can, with certain additional assumptions, provide 
good predictions of the nonlinear results. 
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Figure 14.1: (a) The simulations used to find the point on the manifold je = 1-8, 
K = 15, q/e = 6.3, showing the heat flux vs. q/e at (7^; = 1.8, k = 15). The point 
on the manifold is the point where the heat flux drops to zero, (b-d) Sections 
through the critical manifold with parameters as indicated. Turbulence cannot 
be sustained for k < in (b,c), or for 7^ < 7£;c in (d). The data points were found 
as illustrated in (a), and used to generate the manifold shown in Fig. 114.21 



14.1 Finding the Boundary 

We wish to determine, in a three-dimensional parameter space (7^;, q/e, k), 
the boundary between the regions where turbulence can and cannot be sus- 
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tained nonlinearly. We cover this space using four scans with constant q/e 
(Fig. I14.1l fb)), three scans with constant 7^; (Fig. I14.ir c)) and one scan with con- 
stant K (Fig. 114. ir d)). For each of these cases, we consider multiple values of a 
second parameter and find the value of the third parameter corresponding to the 
zero-turbulence boundary. The boundary is defined as the point where both the 
turbulent heat flux and the turbulent momentum flux vanish. Thus, the location 
of each single point on the boundary is determined using on the order of ten non- 
linear simulations. An example of this procedure is shown in Fig. ll4.lt a). In total, 
we performed more than 1500 simulations to produce the results reported below. 

Because the turbulence that we are considering is subcritical, there is always a 
danger that a simulation might fail to e xhib it a turbulent stationary state because 
of an insufficient initial amplitude llso, 103 1- As we are not here concerned with 
the question of critical initial amplitudes we will consider a given set of param- 
eters to correspond to a turbulent state if such a state can be sustained starting 
with a large enough perturbation. Therefore, all simulations are initialised with 
high-amplitude noise. They are then run to saturation; close to the boundary a 
simulation may need to run for up to t ~ 1000i?/ft/ij to achieve this. 

The critical curves obtained in this manner are plotted in Figs. ll4.lT b-d). These 
curves, which effectively give the critical temperature gradient Kc as a function of 
7£; and q/e, are then used to interpolate a surface, the zero-turbulence manifold, 
plotted in Fig. 114.21 The interpolation is carried out using radial basis functions 
with a linear kernel 196(1 (s^s Chapter ID. 



14.2 The Zero-Turbulence Manifold 

The results of the scan described above are displayed in Figs. I14.ir b-d). These 
three figures show, at fixed values of either 7^;, k or g/e, the threshold in either 
K or q/e below which turbulence cannot be sustained; they are, in effect, sections 
through the zero-turbulence manifold. 

Considering first Fig. 114. It b), we see that, at fixed q/e, the critical gradient 
Kc first rises with 7^;, as the perpendicular flow shear suppresses the ITG-driven 
turbulence, and then falls — in most cases to — as the PVG starts to drive 
turbulence instead. This phenomenon was discussed at length in Chapters |8] and 
[TOl Thus, for every q/e, there is an optimum value of the perpendicular flow shear 
7£; (and hence of the toroidal shear u') for which the critical temperature gradient 
Kc is maximised. We see that reducing q/e increases the maximum that can be 
achieved without igniting turbulence. Fig. I14.ir c) shows that this rule applies for 
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Figure 14.2: The zero-turbulence manifold. Turbulence can be sustained at all 
points outside the manifold (that is, at all points with a higher temperature gra- 
dient and /or higher value of g/e than the nearest point on the manifold). This 
plot is made up from the sections shown in Fig. I14.1l fb-d) (heavy lines) and the 
manifold interpolated from them (thin grey mesh). 



Contours of R/Lt^ 




10 15 20 

u' = dRw / dr / {vthi / R) 

Figure 14.3: Contours of the zero-turbulence manifold plotted against the 
toroidal flow shear u' = dRuj / dr / {vthi/ R) = 7£;/(g/e). The contours indicate the 
value K = Kc below which turbulence is quenched. From top to bottom, the cir- 
cles indicate approximate values of u' and q/e, corresponding to Ref. |l09tl (JET; 
R/L t ^ 8), Ref. [52] (JET ITB; R/ Lj- ~ 17; note large discrepancy, see text) and 
Ref. imh (MAST ITB; R/Lt ~ 10). 
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all considered values of flow sheariH This is to be expected, because lower q/e 
means weaker PVG relative to the perpendicular shear, allowing higher values 
of the perpendicular flow shear to suppress the ITG before the PVG drive takes 
over. 

Lastly, Fig. I14.ir d) shows the threshold in 7^; above which the PVG can drive 
turbulence alone, without the help of the ITG; in other words, even configura- 
tions with a flat temperature profile would be unstable. At very high q/e, already 
a very small flow shear will drive turbulence; as q/e decreases, higher and higher 
values of are required for the PVG turbulence to be sustained. It cannot be 
conclusively determined from this graph whether, as suggested by linear the- 
ory, there is a finite critical value oi q/e below which PVG turbulence cannot be 
sustained, i.e., a nonzero value oi q/e corresponding to 'yEc — ^ 00. However, for 
q/e < 7, the critical 7s is far above what might be expected in an experiment]^ and 
so the 7^; — i- 00 limit is somewhat academic. A definite conclusion we may draw 
is that at experimentally relevant values of shear, pure PVG-driven turbulence 
cannot be sustained for g/e < 7. 

The zero-turbulence manifold interpolated from the numerical data points is 
displayed in Fig. 114.21 The manifold comprises three main features: a "wall" 
where the critical temperature gradient increases dramatically at low q/e; a 
"spur" at low 7^;, jutting out to high q/e (where, as 7^; increases, the ITG-driven 
turbulence is suppressed somewhat before the PVG drive becomes dominant), 
and finally the curve where the manifold intercepts the plane k = 0, whose shape 
is described above. 

14.3 Practical Implications and Comparison with 
Experiment 

In order to illustrate better the implications of our findings for confinement, 
we plot, in Fig. |14.3[ contours of versus q/e and the toroidal flow shear u' = 
dRuj /dr/ {vthi /R)- The basic message is clear: the lower the value of g/e, the higher 
the temperature gradient that can be achieved without igniting turbulence. Once 
we have obtained the lowest possible value of q/e, there is an optimum value 
of u' which will lead to that maximum Kc. We note that the dependence of this 

^ The increase at 7^; = cannot, of course, be due to reduction of the PVG; we assume that 
this occurs because of the simultaneous reduction of the maximum parallel length scale qR in the 
system, leading to weaker ITG turbulence; see [104]. 

^ By order of magnitude, 7^ ^ M/ q, where M is the Mach number of the toroidal flow. Thus, 
values of 7^ much above unity are unlikely to be possible. 
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optimum value of m' on g/e is not as strong as the dependence of the optimum 
value of '^E on q/t (clearly this must be so because u' = {q/t)'^E)- In a device with 
an optimised value of g/e, a near maximum critical temperature gradient would 
be achiev able for u' > 5, shears comparable to those observed in experiment 

lEogl-lmi S. 

While simulation results obtained for Cyclone Base Case parameters are not 
suitable for detailed quantitative comparison with real tokamaks, it is appropri- 
ate to ask whether our results are at all compatible wit h ex perimental evidence. 
For an internal transport barrier (ITB) in MAST, Ref. reports R/Lt ~ 10 

at g/e ~ 4.6 and u' ^ 2.4 ^ This is comparable to the critical values shown in 
Fig. [1431 In JET, Ref. 1I109I] reports R/Lt ~ 8 at q/e ~ 11 and u' ~ 3.6, again rea- 
sonably close to what we would have predicted. However, an ITB in JET studied 



her 



by Ref. |52|1 achieved R/Lt ~ 17atg/e ~ 10 and m' ~ 4.1 — substantially higl 
than our R/Ltc at the same values oi q/t and u' . Note, however, that Ref. |52|l 
reports that the shear was dominated by an enhanced poloidal flow, an effect 
which is not included in our numerical model. 



14.4 Relation to Linear Theory 

Since the mapping of the zero-turbulence manifold using nonlinear simula- 
tions is computationally expensive, we may ask whether linear theory can pre- 
dict marginal stability. The question is also interesting in terms of our theoreti- 
cal understanding of subcritical plasma turbulence. It is clear that in a situation 
where perturbations grow only transiently, existing methods based on looking 
for marginal stability of the fastest growing eigenmode will not be applicable. In 
Chapter [13l we considered these transiently growing modes in a sheared slab, 
and posited a new measure of the vigour of the transient growth: A^max/ the max- 
imal amplification exponent, defined as the number of e-foldings of transient 
growth a perturbation experiences before starting to decay, maximised over all 
wavenumbers. It appears intuitively clear that in order for turbulence to be sus- 
tained, transient perturbations must interact nonlinearly before they start to de- 
cay. We may then assume that a saturated turbulent state will exist if A^max ^ Nc, 
where is some threshold value of order unity. The zero-turbulence manifold 
is then the surface iVmax(7E, k) = Nc- 

^ The ratio of toroidal to poloidal field in MAST can be smaller on the outboard side, so the 
effective value of q/e for locating this case on the zero- turbulence manifold might be smaller than 
quoted. 
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We now test this idea by calculating A'max for linear ITG-PVG-driven transient 



perturbations in a slab, using the code AstroGK [1721] to solve the linearised gy- 
rokinetic equation, as done in Chapter [l3l Fig. I14.4l fa) shows that for each value 
of '^E and a range of q/e, it is possible to choose Nc{'^e) such that the equation 
NmsixilE, f^) = Nc correctly reproduces the critical curve K,c{q/e) obtained as 
a section of the zero- turbulence manifold at that value of 7^. However, Nc does 
have a strong dependence on 7^;, shown in Fig. I14.4r b), ranging from Nc < 0.5 at 
7£;>2toAi'c— ;'Ooas7£;— J-O (the latter is an expected result: at 7^; = 0, there 
is a growing eigenmode, so either A'^max = C)0 or there is no growth at all). It is 
not clear if Nc tends to a finite limit as 7^; — 00, but, similarly to the existence 
of a critical value of g/e as 7^; — )■ 00, this is a somewhat academic question be- 
cause such a limit would be achieved (or not) at 7^; too large to be experimentally 
achievable. 

The practical conclusion of this exercise is that all that appears to be required 
to determine the two-dimensional dependence of Kc on 'Je arid q/e is finding 
Nc{'^e) using a nonlinear scan at a single value of q/e, thus, the number of pa- 
rameters in the nonlinear scan is reduced by one. 




0.5 



1.5 2 2.5 3 3.5 4 

IE 



Figure 14.4: (a) The critical temperature gradient Kc vs. q/e for different values of 
7E, showing both Kc obtained from the interpolated manifold, and Kc such that 



N„ 



Nc, with Nc suitably chosen for each ^e, as shown in (b). 



14.5 Discussion 

In this chapter we have presented two key results. Firstly, and principally, 
we have calculated the shape of the zero-turbulence manifold, the surface that 
divides the regions in the parameter space (7^, q/e, k) where subcritical turbu- 
lence can and cannot be nonlinearly sustained. We have described the shape 

of this manifold and its physical origins, and presented its two implications for 
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confinement in toroidal plasmas: that reducing the ratio q/e, i.e., increasing the 
ratio of the poloidal to the toroidal magnetic field, improves confinement at every 
nonzero value of 7^;, and that at fixed q/e, there is an optimum value of 7^; (that is, 
an optimum value of the toroidal flow shear u' = dRco / dr / (vthi/ R)) at which the 
critical temperature gradient is maximised, in some inst ance s to values compa- 
rable to those observed in internal transport barriers 112 1- How to calculate 



the heat and momentum fluxes that would need to be injected in order for such 
optimal temperature gradients to be achieved was discussed in Ref. | 114| . 

Secondly, we have shown that the zero-turbulence manifold can be param- 
eterised as A^max(7£, 5'/e, ft^) = Nd'yE), where A^max is the maximal amplification 
exponent of linear transient perturbations (calculated from linear theory) and Nc 
must be fit to the data. Thus, using a single scan at constant q/e to determine 
Nc{ie) appears to be sufficient for calculating the full two-parameter dependence 
of the critical temperature gradient. Obviously, the need to fit Nc{'^e) indicates 
a limitation of our current theoretical understanding of the criterion for sustain- 
ing subcritical turbulence in a sheared toroidal plasma. The results reported here 
provide an empirical constraint on future theoretical investigations. 
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Chapter 15 



Summary and Discussion 



We have presented several key results in this thesis. We have used nonlinear 
simulations to show that at zero magnetic shear, a sheared toroidal flow can com- 
pletely suppress turbulence across a range of flow and temperature gradients, 
and can lead to a transport bifurcation from low to high temperature gradients. 
We have identified subcritical PVG-driven turbulence as being a factor limiting 
the temperature gradients that can be achieved in such a transition. We have 
investigated the properties of the transient linear modes which underly this tur- 
bulence, and we have shown that its virulence can be diminished by reducing the 
ratio of the destabilising parallel velocity gradient to the stabilising perpendicu- 
lar gradient (which in a tokamak is equal to q/e). We have used nonlinear sim- 
ulations to calculate, as a function of the parameters (7^;, q/e, R/Lt), the surface 
dividing the regions where subcritical turbulence can and cannot be sustained, 
and from this result we have shown that for low enough values oi q/e, logarith- 
mic temperature gradients R/Lt > 20 can be maintained without generating 
turbulence. 

From these results there naturally arise to further lines of enquiry: practical 
and theoretical. We begin with the practical. From the results given above we 
may infer that tokamak confinement could be greatly increased in the limit of 
low magnetic shear, low q/e and high perpendicular flow shear. However, since 
all three of these quantities must necessarily be a function of the minor radius, 
and since these parameters may only be controlled indirectly in a fusion device, 
it would be necessary to use a 1-dimensional turbulent transport solver such as 
TRINITY [85], with proper models of momentum and heat input, and a more 
realistic magnetic geometry, to demonstrate more conclusively that such a per- 
formance gain could be realised. 

Another avenue for future investigations is determining the dependence of 



the zero turbulence boundary on some of the parameters that were held fixed in 
this work: Tj/Te, magnetic shear, and, more generally, the shape of the flux sur- 
faces, density gradient, inverse aspect ratio e (separately from q), etc. Mapping 
out the dependence just on 7^, q/e and k took approximately 1500 nonlinear sim- 
ulations at a total cost of around 4.5 million core hours. Adding even two or three 
more parameters to the search would take computing requirements beyond the 
limit of resources today, but not of the near future. 

Theoretically, we have demonstrated the existence, at high flow shear, of sub- 
critical turbulence: turbulence which exists in the absence of linear eigenmodes. 
We have described the transiently growing modes which give rise to such tur- 
bulence, and developed a criterion based upon them which had been partially 
successful in determining under which circumstances they may give rise to tur- 
bulence nonlinearly. The theoretical challenge which naturally arises is that of 
developing a full picture of the mechanism of the nonlinear interaction between 
the modes and consequently a robust criterion for when turbulence can and can- 
not be sustained. If such a criterion were to be developed, it would lend the 
ability to model high-flow-shear subcritical turbulence to the quasi linear mixing- 
length-based models which are the workhorse of transport modelling today, as 
well as being a significant advance in our theoretical understanding of subcritical 
turbulence. 

It is clear that there are many questions arising from this work which need to 
be considered. However, these questions aside, we believe that we have provided 
strong evidence that sheared toroidal flows, with a judicious choice of magnetic 
geometry, can be used to achieve high confinement of energy, and hence high 
performance, in fusion devices. 
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Abstract 



The transport of heat that results from turbulence is a major factor lim- 
iting the temperature gradient, and thus the performance, of fusion 
devices. We use nonlinear simulations to show that a toroidal equi- 
librium scale sheared flow can completely suppress the turbulence 
across a wide range of flow gradient and temperature gradient val- 
ues. We demonstrate the existence of a bifurcation across this range 
whereby the plasma may transition from a low flow gradient and tem- 
perature gradient state to a higher flow gradient and temperature gra- 
dient state. We show further that the maximum temperature gradient 
that can be reached by such a transition is limited by the existence, 
at high flow gradient, of subcritical turbulence driven by the parallel 
velocity gradient (PVG). We use linear simulations and analytic cal- 
culations to examine the properties of the transiently growing modes 
which give rise to this subcritical turbulence, and conclude that there 
may be a critical value of the ratio of the PVG to the suppressing per- 
pendicular gradient of the velocity (in a tokamak this ratio is equal to 
qje where q is the magnetic safety factor and e the inverse aspect ra- 
tio) below which the PVG is unable to drive subcritical turbulence. In 
light of this, we use nonlinear simulations to calculate, as a function of 
three parameters (the perpendicular flow shear, qje and the tempera- 
ture gradient), the surface within that parameter space which divides 
the regions where turbulence can and cannot be sustained: the zero- 
turbulence manifold. We are unable to conclude that there is in fact a 
critical value oiqje below which PVG-driven turbulence is eliminated. 
Nevertheless, we demonstrate that at low values of qje, the maximum 
critical temperature gradient that can be reached without generating 
turbulence (and thus, we infer, the maximum temperature gradient 
that could be reached in the transport bifurcation) is dramatically in- 
creased. Thus, we anticipate that a fusion device for which, across a 
significant portion of the minor radius, the magnetic shear is low, the 
ratio qje is low and the toroidal flow shear is strong, will achieve high 
levels of energy confinement and thus high performance. 
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